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ABSTRACT 


This paper contains progress reports of NASA-sponsored studies in the 
areas of space flight theory and guidance theory. The studies are carried on 
by several universities and industrial companies. This progress report covers 
the period from June 15, 1962 to December 20, 1962. The technical supervisor 
of the contracts is W. E. Miner, Deputy Chief of the Future Projects Branch of 
Aeroballistics Division, George C. Marshall Space Flight Center. 
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SUMMARY 

This paper contains progress reports of NASA-sponsored studies in the 
areas of space flight theory and guidance theory. The studies are carried on 
by several universities and industrial companies. This progress report covers 
the period from June 15, 1962 to December 20, 1962. The technical supervisor 
of the contracts is W. E. Miner, Deputy Chief of the Future Projects Branch 
of Aeroballistics Division, George C. Marshall Space Flight Center (MSFC) . 

INTRODUCTION 

Eleven papers are collected in this report. These papers have been written 
by research investigators employed at agencies under contract to MSFC. The 
subject matter lies in the areas of guidance theory and space flight theory. 

This report is the third of the "Progress Reports" and covers the period 
from June 15, 1962 to December 20, 1962. This progress report will herein- 
after be referred to as "the report. " Progress Report No. 1 (2) will be re- 
ferred to as the "first (second) report." Information given in the first and 
second reports is not repeated in this report. 

The agencies contributing and their fields of major interest are: 



2 


Field of Interest Agency 


Calculus of Variations 

Vanderbilt University 
Republic Aviation Corporation 
Grumman Aircraft Engineering Co. 
Hayes International Corporation 
General Electric Company 

Impulse Orbit Transfer 

North American Aviation, Inc. 

Celestial Mechanics 

Republic Aviation Corporation 

Large Computer Exploitation 

Northeast Louisiana College 
Chrysler Corporation 
University of North Carolina 


The objectives of this introduction are (1) to review and summarize each 
agency's contribution to the report and (2) to review the status of each discipline 
and its application to the implementation of the adaptive guidance mode. The 
latter review may be taken as a statement of policy. 

The first paper is written by Dr, M. G. Boyce of Vanderbilt University. 

An application of the theory of the calculus of variations is made to our 
problem. This paper and the note by Dr. R, W. Hunt give enough necessary 
conditions for an optimum trajectory to guarantee sufficiency. The paper repre- 
sents the point of particular interest that the ratios of the Lagrange multipliers 
may be treated as continuous over staging points. The paper is well written 
and most readable. It is recommended for all. 

The second paper is written by Jack Richland of Republic Aviation Corpora- 
tion. A system of equations (deck) for computing "optimum" trajectories is 
presented. Care must be taken to identify all assumptions made by the author. 
The deck is based on fixed end points. The two control variables are thrust 
direction and thrust magnitude. The mass at injection (final cutoff) is maxi- 
mized. The deck permits the design of trajectories for vehicles with restartable 
engines. It may be used for fuel loadings in given stages. A differential cor- 
rection method is used to establish initial values of the Lagrange multipliers. 

The method should work well with good initial guesses. The fixed end-point 
assumption is not overly restrictive. How r ever, many problems require func- 
tional relationships for injection conditions. In this case the fixed end point 
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solution requires added optimizations. There are computational advantages to be 
gained by not using the assumption. Because transversality conditions may be de- 
rived from injection conditions, these conditions assure an optimum final answei. 
The differential correction method described in the report may be used with 
functionally defined injection conditions. The transversality conditions, in- 
cluding those for variation in burning time and mass, should be used to adjust 
[3> (tp) ] as is discussed on page 4? . 

The fourth and fifth papers present methods for computing low-thrust tra- 
jectories in the neighborhood of a large attracting body. The fourth paper, by 
Gordon Pinkham of Grumman Aircraft Engineering Corporation, presents formu- 
lations for direct methods of the calculus of variations. The fifth paper, by 
Harry Passmore of Hayes International Corporation, presents formulations for 
classical methods of the calculus of variations. The parameters used in both 
papers are basically the same. The differences are in the way the parameters 
are handled. Pinkham develops the differential equations for the parameters. 

The Euler- Lagrange equations for the problem are also found. These are 
given as equations 1, 2, and 3 along with equations 5. Equations 5 and 6 give 
the basic equations for the gradient procedure. This procedure has been 
checked for feasibility. The potential is considered good. The Passmore paper 
takes as a base the planetary theory of Lagrange. The vehicle is taken as one 
planet. The thrust action is taken as representing a second planet. The coupling 
effect is the perturbation for each planet. The perturbation differential equa- 
tions are then developed in standard form. These will be integrated and the 
results will be analyzed before further development is attempted. 

The next paper, by Carlos R. Cavoti of General Electric Company presents 
a problem with a restricted control variable. The assumptions are: (1) The 
rocket is always in vertical climb, (2) the thrust magnitude, bounded from above 
and below, is the control variable, and (3) the non-potential forces considered 
are those of thrust and drag. Drag is considered to be a function of velocity and 
altitude. The paper develops methods of finding corners, i. e. , the points at 
which the thrust gradient varies. A problem is used to illustrate the procedure. 

The paper by Dr. D. F. Bender of North American Aviation covers a 
special field in optimization. Here only impulsive forces are considered. This 
study follows extensive work on orbital transfer. It may also be considered as an 
early investigation on rendezvous. This work and the techniques developed are 
needed for early trajectory design. 

In the field of the calculus of variation, the studies in theory for high thrust tra- 
jectory problems are sufficient for our present needs. Extended effort needs to 
be expended in three other areas, however. The first of these areas is that of 
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low thrust trajectory calculations. The second area is in the development of 
specific decks for special problems. The third area is in developing methods for 
expanding the A values at t Q (t). These expansions will be in terms of the stated 
variables at t Q (t) and the desired end conditions (including motor characteristics). 

Two papers in this report cover the area of optimum low thrust trajectories. 
Both papers present feasible methods. Grumman Aircraft Engineering Corpora- 
tion and Hayes International Corporation will continue to work in these fields . 

Auburn University has been developing decks for specific problems. 

They will continue to do this type of work. Such work requires an eye for 
mathematics, physics, and hardware. Other contractors may be added to assist 
in these studies as the need arises. 

The third ai'ea is the -least explored. There are two potential methods 
of making the expansions for Aj ( t 0 ) . The first is numerical; the second is 
analytical. 

An outline of an approach for numerical expansion is given below. Refer- 
ence is made to the comments on the Republic Aviation paper. The first partial 
derivatives have been developed in a manner similar to the one described. The 
transversality conditions are related to the Euler- Lagrange equations. The 
Jacobi equations are of one higher order. Thei'e are conditions (say, second 
difference transversality conditions) related to the Jacobi equations which, 
when combined with the third partial relations, (say, second Jacobi equations) 
give the second partial derivatives. Higher partial derivatives may be developed 
in a similar manner. Dr. R. W. Hunt has developed the second transversality 
conditions. This work should proceed with accuracy and speed. The differential 
generator developed by the University of North Carolina will be used to develop 
the partials. This will be done at Marshall Space Flight Center. The contributors 
will develop the operations needed for use of these differentials only. The steps 
needed in the development and the contractors involved are given below: 
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Step 

Development 

Contractor 

1 

1st transversality condition 

completed 

2 

1st Jacobi equations 

completed 

3 

2nd transversality condition 

Dr. Hunt (SIU) 

4 

2nd Jacobi equations 

Dr. Boyce (VAN) 

5 

3rd transversality condition 

Dr. Boyce (VAN) 

6 

3 I'd Jacobi equations 

unassigned 

7 

4th transversality condition 

unassigned 


The analytic approach will not be detailed. 

The end conditions are 


defined by the mission criteria formulation and the transversality conditions. 

We will assume a formal expansion for the solution of the differential equations; 
Chapter II of reference 2 gives such an expansion. The differential equations 
solutions will be substituted into the end conditions. The result will be infinite 
series in time at cut-off. The coefficients of the expansion will be functions 
of initial and end conditions. These expansions will be inverted for the initial 
A's and time at cut-off. This inversion may be attempted by formal means. 

Only the University of North Carolina and Marshall Space Flight Center have 
considered this problem to date. Other contractors may be assigned special 
parts of this problem. 

The impulse optimization by North American Aviation will be continued 
until this field is thoroughly covered. 

There is only one paper on celestial mechanics. It is by Dr. Mary Payne 
of Republic Aviation Corporation. In this paper the origin of the coordinate 
system is accelerated and the coordinate system is rotating. The coordinates 
are the initial values in the solution of Euler's problem of two fixed centers. 

The origin is selected in such a way as to minimize the effects of non-integrable 
terms in the perturbation equations. Work will continue on this problem to 
determine the value of the approach. 

There are three papers in the field of large computer exploitation. The 
first of these is by the Northeast Louisiana State College Department of Mathe- 
matics group. This paper presents a recursion process which is to be used in 
least squares cux*ve fitting. In the process the points taken from the trajectory 
are designated by a parameter where /3j is the ordered m-tuple , 
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1^8’ tii, • • • ^im ]. There seems no need for this restriction. The vector 
[ti 0 , ty, . . . ti m ] may be any m functions such that x may be considered a 
function of these functions. This is implied by the last paragraph. 

We take as given values ty (i = 1 . . . n, the points; j = 1 . . . n (functions) ) 
and xp Define 

1) gj= tfj i = 1 . . .m 


2) s' y = g y - (g y • e 0 ) e 0 - ( g^ • e t ) e t - .. 
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X 


Si 


7 



e ) e 

y-i y-i 


Let y successively take on values y = 0, 1, . . . , m. and compute g^, || g^ 

and e* where g f = g a 

y °0 u 

3) Define x = [XivvX n l 


4) A. = x • e. j = 0, . . . , m. 

J J 

Mr. Vance of Chrysler presented a method for solving the above problem in the 
second progress report. It differs in outlook and procedure from the above. 


The present report by Mr. Vance of Chrysler purports to present M one way 
of solving the problem of data point selection for the generation of a least 
squares approximation of a multivariate function by a linear combination of 
polynomials which are orthononnal over a region. " It is needless to say this 
problem needs our attention. The work developes only the concepts. Much added 
work is required before the method is usable. Time has not permitted the search 
of references. As a result certain questions remain. Some of these are: (1) 
what are the details of the derivations of equations on page 1915(2) what modifi- 
cations are needed to replace Wj with density considerations? Mr. Vance will 
continue to work in this area. The paper shows potential for success in developing 
a means of selecting points for least squares curve fitting. 
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The last paper of this report is by Shigemichi Suzuki and Sylvia M. 

Hubbard of the University of North Carolina. Linear programming techniques 
are used to fit known functions. The linear programming techniques should better 
control the errors at injection. The work presented is to be implemented at 
Marshall Space Flight Center in the near future. 

There are several unanswered questions in the multivariant functional 
model development. Four of these questions are listed. 

1) What terms should be selected for the polynomial? 

2) What criteria should be used in developing the polynomial? 

3) What points should be selected for fitting the polynomial? 

4) What is the best functional form for fitting the data? 

Suzuki and Hubbard’s contribution will let us start to look at questions 
two and four. Vance's work attacks the third question. No theoretical work 
has been done on the first question. 

The accomplishments in this field have been large. The least squares 
approach has been used. Polynomials have been used as the functional form. 

This gives an engineering answer to questions two and four. Empirical methods 
have been used to answer questions one and three. The results are sufficient 
for the needs of today. 

A complete steering angle program and a time-to-fly program have been 
developed here in the past for a specific vehicle. The time required for the 
development of each step is recorded below. 

a. One week to gather data if data is readily available. 

b. Three days to establish the performance for the 1st stage. 

c. One day to obtain the volume of 1st stage trajectories. 

d. One day to obtain the volume of second stage trajectories and invert 
the matrix, (f" hr for calculus of variations plus <r hr for transferring 
data plus one hour for matrix inversion for a total of 2 hrs machine 
time) . 
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e. One day to compute coefficients ( 15 minutes of machine time). 

f. Two days to record and check coefficients. 

g. One day to check selected x curve by simulation runs ( l|- hour machine 
time) . 

The above were actual times required. The large time for small computer 
runs was caused by data handling and delays between computer runs. It will be 
noted that two days were lost in recording data. Two days were also lost by 
human error. Automation of the process is underway. There is no known way 
to reduce the data gathering time. Complete automation should make the re- 
mainder of the problem an over-night job. It will also cut the cost of human 
error. Automation of the process will make procedures static. It is costly 
and it should be done only when the state of the art permits. Whether the 
state of the art permits such automation now is an open question. 


REFERENCES 

1. Hunt, Robert W. , "Utilization of the Accessory Minimum 
Problem in Trajectory Analysis," MTP-AERO-62-74, 

October 2, 1962. 

2. Moulton, F. R. , "Differential Equations," Dover Publishing, 
Inc . , 1930. 
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AN APPLICATION OF CALCULUS OF VARIATIONS 
TO THE OPTIMIZATION OF MULTISTAGE TRAJECTORIES 


By 

M. G. Boyce 


SUMMARY 

A procedure is developed for determining a fuel minimizing trajec- 
tory for a multistage rocket in three dimensional space. In each stage 
the fuel burning rate and magnitude of thrust are assumed constant. The 
motion is subject to the inverse square gravity law but with negligible 
atmospheric resistance. The Euler -Lagrange equations determine minimizing 
trajectories in a given stage. Transversality conditions are then invoked 
to extend a minimizing path across the boundary to the next stage. The 
existence of minimizing trajectories is assumed, sufficient conditions 
not being investigated in this paper. 

In Appendix I a simple form of Zermelo's navigation problem, extended 
to several stages, is solved to illustrate some aspects of multistage 
problems . 

Appendix II gives a summary of necessary conditions for calculus of 
variations problems of the Mayer form involving control variables. 
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I . INTRODUCTION 

The problem is to determine the fuel minimizing trajectory of a 
rocket whose flight consists of several stages caused by engine shutoffs 
at specified times. Initial position and velocity are assumed given and 
target conditions specified. In each stage the analytic formulation is 
similar to that of Cox and Shaw (Ref. l), and we make their basic assump- 
tions that the earth can be considered spherical, the inverse square 
gravity law holds, the only forces acting on the rocket are thrust and 
gravity, the direction of thrust is the axial direction of the rocket, 
rotation effects can be ignored, in each stage the magnitude of thrust 
and the fuel burning rate are constant, and the center of mass of the 
rocket is fixed with respect to the rocket. 

The general procedure is roughly as follows. Using the fixed initial 
conditions for the first stage, determine as solutions of the Euler-Lagrange 
equations the family of minimizing trajectories satisfying those conditions. 
The given time t x for the end of the first stage will fix on each mini- 
mizing trajectory a definite point. The totality of these points will 
constitute a subspace S , which will be the locus of initial points for 
the second stage. New values of mass, thrust, and fuel burning rate de- 
termine new Euler-Lagrange equations. Minimizing trajectories must satisfy 
these new equations in this stage and also must satisfy transversality 
conditions for initial points in subspace . Through each point of S.^ 

these conditions determine a unique trajectory, and on each of these tra- 
jectories the given time t 2 for the end of the second stage will fix a 
definite point. The totality of these points will be subspace S 2 , which 
in turn will be the locus of initial points for the third stage, and trans- 
versality conditions will again determine a family of minimizing trajec- 
tories, one issuing from each point of S g . This procedure is repeated 
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until in the final stage the mission objectives will impose criteria for 
selecting a pieced trajectory satisfying the given initial conditions and 
extending through the several stages. Closed form solutions are not 
attainable in most cases. However, it would seem possible to extend the 
single stage adaptive guidance mode computational procedures through 
several successive stages. 


II. FORMULATION OF THE PROBLEM 


A plumbline coordinate system is used (Ref. 1, p.108; Ref. 2, p.ll), 
with the center of mass of the rocket designated by x = (x ,x ,x ) and 
its velocity by u = (u^ > u 2 > u 3 ) • The time t is taken as independent 
variable, and u = dx/dt.. The thrust vector F = (0,F,0), having its 
magnitude F constant for each stage, is assumed to be directed along 
the axis of the rocket. The orientation of the rocket axis relative to 
the plumbline system is designated by ^ , where % ± ,% 2 ,% 3 

are the pitch, roll, and yaw angles, respectively. 

If denotes the gravitational acceleration and [a] the matrix 

for transformation of vectors from the missle to the plumbline coordinate 
system, then Newton's second law gives as equations of motion of the rocket 


u = m" 1 F [aJ + x , x = u . (l) 

o 

In terms of pitch, roll, and yaw, the matrix A has the following form 
(Ref. 1, p.108; Ref. 2, p. 26 ): 


A 


CPCR SPCR SR 

-SPCY - CPSRSY CPCY - SPSRSY CRSY 

SPSY - CPSRCY -CPSY - SPSRCY CRCY 


CP = COS “/q 
SP = Sin % 1 
etc . 


Since roll effects are to be ignored, the roll will be assumed 

identically zero. Hence CR = 1, SR = 0, and the variable % 2 may be 
dropped. Since fuel consumption is monotonically increasing with time, 
minimization of time of flight is equivalent to minimizing fuel consump- 
tion. It is more convenient to treat the problem from the minimum time 
standpoint . 
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In the terminology of the genera] theory of Appendix II we now have 
state variables u ,u 2 ,u 3 , x 1 > x 2 > x 3 > control variables and % 3 , and 
independent variable t. The function to be minimized, the function 
h(b) in the appendix, is simply the final time t f . Hence t f is one 
of the parameters in b ; other parameters may occur in the initial and 
end conditions and in stage boundary conditions. The mass m is assumed 
a known function of t in each stage so is not included in the state 
variables . 

Thus the problem is to find in a class of admissible sets of functions 
u(t), x(t), %(t) and parameters b a set that will satisfy the differen- 
tial equations (l) and the given end conditions and that will minimize 
the final time t . 


III. FIRST STAGE 


Let the time interval for the first stage be t ^ t t^ and the 


initial conditions, u(t ) = 


— o 


u 


X 


— o — o 


(t ) = x . On putting % s 0 in A 


and using jur" 3 x for x , where (u is the gravitational constant times 
the mass of the earth, we get equations (l) in the form 
u, = -Fm _1 SPCY - ur" 3 Xi 


u = Fra -1 CPCY - ur"°x„ 
2 .2 

u„ = Fm" 1 SY - ur" 3 x„ 

X = U 
1 1 

*2 = U 2 


X 3 = U 3 


( 2 ) 


In order to apply the necessary conditions of Appendix II, we now 
define a generalized Hamiltonian 

H = L 1 (-Fm _1 SPCY - |ar" 3 x i ) + L 2 (Pm“ 1 CPCY - |ur~ 3 x 2 ) 

+ L g (Fm -1 SY - ^r" 3 x 3 ) + + L 5 u £ + L 6 u 3 . 

By condition I, Appendix II, the Euler-Lagrange equations are 


U - 


h >h mU L > \ - -“u. ’ U - -t. ’ y = °’ 1 ° 1 ’ 2 ' 3 ’ 

1 1+3 1 1 J J = 1,3. 
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These formulas give the six equations (2) plus the following eight: 



3R- 5 x i (L i x i + L s x 2 + L 3 x 3 ) 
3 F r- 5 x 2 (L i x 1 + L 2 x 2 + L 3 x 3 ) 
3 F r- 5 x 3 (L iXi + L 2 x 2 + L 3 x 3 ) 


(3) 


L CPCY + L SPCY = 0 
1 2 

L SPSY - L CPSY + L OY = 0 

12 3 


(4) 


Assuming CY ^ 0 and letting D 2 = L 2 + L 2 , E 2 = L 2 + L 2 + L 


i " 2 

we get from equations (4) that 

tan q = -L i /L 2 , SP = -L x /D , CP = L £ /D , 

tan i 3 = L 3 /D , SY = L 3 /E , CY = D/E , 


1 2 
D > 0, E > 0 


(5) 


the choice of signs in SP,CP,SY,CY being a consequence of the Weierstrass 
and Clebsch conditions, as will be shown in the next section. From ( 5 ) it 
follows that the thrust vector in the plumbline system can be expressed as 


F [a] = F ( -SPCY , CPCY , SY ) = F(L /E,L 2 /E,L 3 /E) . 

Equations ( 5 ) ma y be used to eliminate the control variables from 
equations ( 2 ), thus giving, together with equations ( 3 ), a system of 12 
differential equations of the first order in 12 dependent variables. 
This system may be written as six equations of second order, which in 
vector notation are 


x = FE/mE - jix/r 3 , 

E = -uE/r 3 + 3u(x-E)E/r 5 ; 
where E denotes the vector (L ,L ,L ) . 

“ 1 7 2 7 3 7 


(6) 
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Although the result is not utilized in this paper, it is of interest 
to note that three first integrals of the system (6) can be readily ob- 
tained by the following device. Cross multiply the first of equations 
(6) by E and the second by x and add the resulting equations to get 

E)(x + x/ E= °. ^ 

This now yields 


EXi + x/ E = M, ^ 

where M is a constant vector, since the derivative with respect to t 
of the left member of (8) is the left member of (7). 

The equations (2) and ( 3 ), after elimination of the control variables, 
or, equivalently, system (6), will have a six -parameter family of solu- 
tions satisfying the given initial conditions u(t Q ) = H o ,x(t Q ) = - Q • 
However, since the equations are homogeneous in the L's, if u(t),x(t), 
L(t) is a solution, then so is u(t ) ,x(t ) ,cL (t ) for any non-zero con- 
stant c . Thus, if initial values of the L's are taken as parameters, 

only their ratios are significant in determining u(t),x(t) . Hence the 
value of one L may be fixed, or some function of the L's may be as- 
signed a value at t = t Q , say LptJ + L|(t Q ) + L^(t Q ) = 1 . Thus there 

is a five -parameter family of trajectories satisfying the Euler -Lagrange 
equations and having the given initial values. If b_^,...,b^ denote the 
parameters , the equations of the family may he written 


u = u(t,b x ,b 2 ,b 3 ,b 4 ,b 5 ) , (9) 

X = X (t ,b ,b^ ^3 ^4 5 ^ " 


Each of these curves is the path of least time from the initial 
point to any other point on it, assuming that a minimum exists and that 
only one of the curves joins the two points. (The geometrical terminology 
refers to the seven dimensional space t,u,x and not to three dimensional 
physical space.) Putting t = t x gives a point on each curve, and the 
totality of such points constitutes a sub space S 1 . If S 1 is con- 
sidered as a given locus of variable end-points for the first stage, then, 
since t has constant value p on S 1 , each trajectory is a time mini- 
mizing trajectory from the initial point to S 1 , and hence must satisfy 
the transversality conditions at . This property will be utilized in 

Section VI. 
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IV. THE WEIERSTRASS AND CLEBSCH CONDITIONS 

We now show that, with the choice of signs adopted in ( 5 ), the neces- 
sary conditions II and III of Appendix II are satisfied by solutions of 

equations ( 2 ), ( 3 ), (4). For the Weierstrass test let circumflexes denote 
arbitrary values of the control variables. Then 

H(t,u,x,2^,L) “ H(t,u,x,7C,L) 

= Em" 1 (-L i SPCY + L 2 CPCY + LgSY + I^SPCY - L 2 CPCY - L 3 SY) 

= Em -1 (E + L SPOT - L 2 CPCY - L 3 SY) > 0 , 

as is implied by the general inequality 

(a 2 + b 2 + c 2 ) 1/2 > (d sin A + b cos A) cos B + c sin B , 

which holds for all real values of a, b, c. A, B. 

For the Clebsch test , the matrix of the quadratic form involved is 

LjSPCY - L 2 CPCY L 1 CPSY + L 2 SPSY ] 

I^CPSY + L 2 SPSY L 1 SPCY - L 2 CPCY - L 3 SY I . 

By virtue of equations ( 5 ) this becomes 

-D 2 /E 0 

0 -E 

l J 

which implies that the quadratic form is negative definite. 

There are in all four sets of values of SP, CP, SY, CY in terms of 
the L's that will satisfy equations (4). Two of them reverse the in- 
equality signs in conditions II and III, but there is one other set besides 

that given in ( 5 ) that satisfies conditions II and III. It can be got from 

( 5 ) by replacing D by -D . This amounts to changing to ^ + it 

and X 3 bo it - ^3 > an( ^ is found that this actually produces the same 
direction of thrust as before. 
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V. SECOND AND SUBSEQUENT STAGES 

For the second stage the range of t is t^ < t < . The initial 

point is required to be in S^ , the equations of which are ootained by 
putting t = t in ( 9 ): 

u = u(t ,b ,b ,b 3 ,b 4 ,b 5 ) "\ 

[ - X x (b) , (10) 

X = x(l ^ ^ ’^2 ^3 ^ «/ 

the six functions in the right members being denoted by X^) to con- 
form with the notation in Appendix II. The function T^b) i s the 
constant t 

1 

The differential equations of motion are of the same form as for the 
first stage , although F and m have different values. To allow for 
possible discontinuities in the L's , we denote their right hand limits 
at t by L(t i +) . There are five transversality conditions (Condi- 
tion I, Appendix II ) which must be satisfied at t = t^ : 

L(t +) • X , = 0 , k = 1 , 2 , 3, 4, 5- (H) 

— x 1 —ib. 

k 

Since these equations are homogeneous in the L's , and so are the equa- 
tions analogous to ( 2 ), ( 3 ), and (4), it follows that for the determina- 
tion of u(t) and x(t) again only the ratios of the L's are signifi- 
cant. Thus again there will be an eleven parameter family of minimizing 
trajectories. When values are given to the b's to fix a point in S 1 , 
there will be six values u(t x ) , x(t 1 ) and five transversality conditions 
to determine the eleven parameters. This in general will fix a unique 
minimizing trajectory issuing from each point of S^ . Let the equations 
of these trajectories be expressed by the same equations ( 9 ) as for the 
first stage except that now the range for t is from t 1 to t 2 
Putting t = t will determine a definite point on each trajectory, and 

the locus of these points will be a subspace S g with equations 
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u = u(t .b .b .b .b .b_) 

— — 2 * 1 * 2 * 3 * 4 * 5 ' 

x = x(t ,b ,b .b .b .b ) 

- - 2 * 1* 2 9 3 9 4 9 s' 


= x 2 (b) # 


Note that again the transversality and other conditions involving the end 
point need not be used to determine the five parameter family of trajec- 
tories but only the conditions at the initial point. 

For subsequent stages the procedure is like that for the second 
stage. The initial point for the third stage would be restricted to 
subspace and transversality conditions involving X 2 (b) and 

L(t +) would be used. 

The computational procedure given by Cox and Shaw (Ref. 1, pp 118) 
could be used in the first stage. Modifications would be needed in the 
other stages to approximate the partial derivatives of the X(b) functions 
and to solve the transversality equations. 

In the final stage the mission objectives must be fulfilled at the 
end point. Since there is little hope for closed form solutions, the 
proposed procedure is to estimate initial conditions and use them to ex- 
tend a solution by approximate integration methods through the several 
stages. If the objectives are not attained, make new estimates of initial 
conditions and new computations of a minimizing trajectory, continuing 
thus until a trajectory is obtained that achieves the desired objectives 
with sufficient accuracy. 


VI. CONTINUITY PROPERTIES OF THE LAGRANGE MULTIPLIERS 

In each stage the trajectories which are without corners and which 
satisfy the Euler-Lagrange equations will have Lagrange multipliers that 
are continuous and differentiable (Ref. }, pp 202-204; Ref. 6, pp 12). 
However, on passing from one stage to the next, there are discontinuities 
in the functions defining u . From equations (2) it follows that there 
will be corners for the functions u , and hence discontinuities might be 
expected in the L's . But the functions defining x and L are con- 
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tinuous in t, u, x and have continuous partial derivatives. Thus con- 
tinuous solutions for the L's can be obtained by taking u continuous 
across boundaries, provided the transversality conditions can be satisfied. 

In obtaining the family of solutions of the Euler -Lagrange equations 
in each stage the homogeneity of the equations in the L's was utilized 
to decrease the number of parameters by one, say by assigning an initial 
value to one of the L's . As remarked at the end of Section III, the 
five transversality conditions for parameters ‘ b 1 >“*> b 5 , namely, 

L(t^-) ‘ ( b ) = 0 > k = 1,2, 3, 4, 5, 

are satisfied on S i . These conditions are the same as conditions (ll) 
in L(t +) which hold for S ± as locus of initial points in stage two. 
Hence and ^ (t+) , . . . ^(t^) are proportional. 

By assigning equal values to one pair from the two sets, all can be made 
continuous at t . 

The transversality condition involving the final time as parameter 
in each stage is not homogeneous in the L's because of the term h^ 

k 

This condition would make the set of L's unique and not necessarily 

continuous across the boundary; however, it is not essential to use this 

condition for the determination of the trajectory equations. Hence it is 

possible to obtain Lagrange multipliers that are continuous through the 

several stages and to use their ratios at the initial point t — t^ as 

parameters b , . . . ,b for a five parameter family extending through all 
L 1 5 


the stages. 
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APPENDIX I 

A MULTISTAGE NAVIGATION PROBLEM 

A simple form of Zermelo's navigation problem (Ref. 4), extended to 
multiple stages, serves to illustrate some features of trajectory problems. 
Zermelo stated his problem for air flight in a plane, but we follow Cicala's 
formulation (Ref. 5 , pp 19) and consider a motor boat on a plane water 
surface. A rectangular coordinate system is associated with the plane 
surface, and the boat is considered a point (x,y) . The water current 
is assumed to have known velocity components u and v as functions of 

x and y and the time t . Let the velocity vector of the boat relative 

to the water make an angle 0 with the positive x-axis and assume that 
the magnitude of the velocity vector is a known constant in each stage. 

The path of the boat is determined by the control variable 0 , and the 
problem is to find 0 as a function of t so as to minimize the time 

for the boat to go from the origin to a specified point (x f ,yj that 

is assumed remote enough to require three stages. In order to get a prob- 

lem that vi ill have an easily obtained closed form solution, we take the 
water velocity components to be constants and choose the coordinate system 
so that u = 0 , v = a . 

The problem then is to find functions x(t), y(t), 0(t) such that 

x = v cos 0 , y = a + v sin 0 ; (l) 

v = v i for 0 < t < t i ; v = v- 2 for ^ < t < t g ; v = v 3 for t < t; 

x(0) = y ( 0) = 0 ; x(t f ) = x f , y(t f ) = y f ; 
and such that t^ is a minimum. 

First Stage . 

As in Appendix II, define the generalised Hamiltonian 
H = L v cos 0 -i- L (a + v sin ©) . 

11 2 i ' 
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From this H the Euler-Lagrange equations are Found to be 


x = v cos 6 , y = a + v sin 0 , 

. 1 . ' ( 2 ) 

L = 0 , L = 0 . -L sin 0 + L cos 0=0. 

1 - n 7 1 2 

Hence and .are constants , say = L^, L p . It then 

follows that e is constant, and integration of the first two of the 
above equations gives 

x = (v eos 0)t , y = (a i v i sin 6)t , (o) 

on using initial conditions x = y = 0 when t = 0 . Thus paths of mini- 

mum time are straight lines. 

If our problem were a one -stage problem with end point (x ,y ) and 

time t to be a minimum, we would have for the determination of 9, 

L and L the following equations 
11 21 

x = (v cos 0)t , y=(a+vsin0)t, (^ ) 

l i li i i 

-L sin 6 + L cos 6 = 0 , (5) 

11 21 

plus the transversality condition 

L v cos ® + L 2i^ a + v i sin - 1 ■ (6) 

Equation (6) is found from the transversality equation in Appendix II by 
putting 


* = 1 > h - V T 1 - °’ T 2 - V X 11 = °’ X 21 = °’ X 12 = V X 22 ’ y i> h " V 

Equations (b) determine 0 and t , while (5) and (6) give unique multipliers 

L = cos 0/(v + a sin 0) , L = sin 0/ (v + a sin 0) . (7) 

11 1 2 1 i 


Now if we consider ^ ,y ) variable and inquire as to the locus of 

such points each of which is reached in a minimum time equal to t , we 

get from (4) with 0 variable that the locus of (x ,y ) is the circle 

with center (Chat ) and radius v t 

i li 



Second Stage . 

The locus of initial points for the second stage is the circle men- 
tioned in the preceding sentence. We write it as 

x = (v cos a)t , y = (a + v sin cl ) t (8) 

11 7 i ; i l 

with the parameter OL replacing the 0 of equations ( b ) since we shall 
continue to use 0 as the control variable. The differential equations 
of constraint for this stage are the same as for the first stage except 
that v 2 replaces v . 

The Euler -Lagrange equations are as before, with v^ replacing v , 

and hence L and L are constant, say L = L , L = L . It follows 
1 2 7 " 1 12 7 2 22 

that 0 is constant. 

If the end point for the second stage is considered fixed at ^x 2 ,y 2 ) , 
then transversality conditions for parameters OL and t are 

L v t sin OL - L v t cos OL = 0 , 

12 1 1 22 1 1 

, , (9) 

L v cos 0 + L„(a + v sin 6) = 1 . 

12 2 22 2 

The first of these equations, together with the last of the Euler -Lagrange 
equations, implies that 0 = OL . Then, from the pair of equations (9)> it 
follows that 

L =cos 0/(v 2 + a sin 0), L^ =sin 0/ (v^ + a sin 0) . (10) 

Thus L and L are not equal to L and L , indicating di scon- 
12 22 1 11 21 0 

tinuities in the multipliers at stage boundaries. However, the control 

variable 0 is continuous, being in fact the same constant in the two 

stages . 

On integrating the Euler -Lagrange equations for x and y and using 
(8) as initial conditions, one finds that 

X — ( v cos 0)t + (v - v )t cos 0 , 

y = (a + v 2 sin 0)t + (v jL - v 2 )t sin 6 . 


( 11 ) 
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For each constant 0 , the path is a straight line. 

Wow consider the locus of end points (x 2 ,y 2 ) that will each be 

reached in minimum time t £ . Fixing t = t 2 in (ll) and considering 
0 variable shows the locus to be the circle with center (0,at 2 ) and 
radius v t.^ + v 2 ^ 2 " ‘ 

Third Stage . 

With the circle of the preceding sentence as locus of initial points , 
the end point is required to be (x^,y^) and time t^ is to be a minimum. 

In the same way as before the path is shown to be a straight line with the 
control variable constant and equal to its value in the preceding stages . 

The new equations for x and y are 

x = (v cos 9 ) t + f (v - v )t + (v - v )t 1 cos 9 , 

3 1121 2 j (12) 

y = (a + v 3 sin 0)t + [ (v i - V 2 H 1 + ( v 2 ” V 3 ^ 2 ] s ^" n ® ‘ 

By putting the given values x^,y^ in equations (12), one can solve for 
the minimum time t = t~ and for the constant control angle 0 . Then 
equations (ll) with t = t 2 , x = x g , y = y 2 and equations (8) determine 

the corner points (x^y^ and (x 2 ,y 2 ) . 

Conclusions . 

This problem illustrates the extension of a trajectory across stage 
boundaries where the differential equations of constraint are discontinuous. 
The effect of the homogeneity in the Lagrange multipliers is similar to that 
in the more general problem. 

The unique Lagrange multipliers that satisfy the Euler -Lagrange equa- 
tions and the transversal ity conditions of I, Appendix II, are discon- 
tinuous at stage boundaries. However, the ratio L^L^ = tan 9 is the 

same for each stage. The equations containing L's are homogeneous in 
the L's , except that the transversality condition computed for the final 
time as parameter in each stage is not homogeneous. But this transversality 
condition is not needed to determine the family of minimizing trajectories 


which satisfy initial conditions in each stage. That is, in order to ob- 
tain a pieced trajectory extending through the several stages, only the 
ratio of the L's is needed, and, since the ratio is preserved, the L's 
may be chosen continuous. 

The geometrical interpretation of this problem is of interest, so a 
diagram is shown below. The computations were made for 

a = 1 , t i = 2, t 2 = 5, v 1 = 3, v 2 = 1 , v 3 = z, x f = 15, y f = 10, 
and yielded the results 

t f = 8 . 09 , 9 = 7°lV, (x x ,y^ ) = ( 5 . 93 , 2.76), (x 2 ,y 2 ) = (8.89, 6.13). 

\ 



( 10 , 15 ) 




APPENDIX II 


NECESSARY CONDITIONS FOR MAYER PROBLEMS 
WHICH CONTAIN CONTROL VARIABLES 


This appendix lists the principal classical necessary conditions, 
modified to allow for control (or undifferentiated) variables in the 
constraints. For proofs and fuller discussion see References 6, 7, 8* 


Notation 

x = (x 1? . . . ,x n ) 

y = 

b = (b 1 , . . • >b_^) 

R-Si = (Li—'V 
= (x 21 ,...,x 2n ) 


g “ (g x ? * * * 

L = (L^ , • • • jL^) 


H = L • jg 


state variables , functions of independent variable t 

control variables , functions of t 

parameters occurring in end conditions 

functions of b defining first end point 

functions of b defining second end point 

functions of (t,x,y), defining derivative constraints 

Lagrange multipliers, functions of t 

generalized Hamiltonian function 


h{b) function to be minimized 

Variables occurring as subscripts will denote partial derivatives, 
and a superimposed dot will indicate differentiation with respect to t. 

A set t,x,^,b will be called admissible if it belongs to a given open 
set R , and a set x{t) ,y(t) ,b will be an admissible arc if its elements 
are all admissible and if x(t) is continuous and x(t),y(t) are piece- 
wise continuous. The functions occurring in T ; X ; _g, and h are assumed 
to have continuous partial derivatives of at least the second order. 
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Statement of Problem 

In a given class of admissible functions and parameters x(t),£(t),b 
it is required to find a set which satisfies the differential equations 
and end conditions 

x = g(t ,x,^) , t i ^ t ^ t 2 

t x = T x (b) , t 2 = T 2 (b) , x(t x ) = X x (b) , x (t 2 ) = Xg(b) 

and which minimizes the given function h(b) . 

Let C be an admissible arc x(t), y( t), b which is a solution of 
the problem. Also let C be assumed normal (Ref. 6 , pp 15 ) and to have 
’x(t) and £(t) continuous. Then C must satisfy the following four 
conditions . 


I. First Necessary Condition . For every minimizing arc C there 
exist unique multipliers L_.(t), having continuous first derivatives, such 
that the equations (Euler-Lagrange) 


Hj-^ , L. — H , H — 0^ i — l,...,n, j — l,...,m, 

i i 

hold along C . Also the end values of C satisfy the transversal ity 
conditions 


V b, - h. ■ h K - VW + L • 

k k k 

where subscripts 1 and 2 on H and 
t = t and t = t } respectively. 

As a consequence of the above Euler 
that also along a minimizing arc C 


+ ix b = °* k = l,---,r, 
k k 

L indicate evaluation for 


-Lagrange equations it follows 


dH/dt = H , 

and hence that, if H does not involve t explicitly, then H is con- 
stant along C . 


II. Weierst.rass Condition . Along a minimizing arc C the inequality 
H(t,x,Y,L) < H(t,x,£,L) 

must hold for every admissible element (t,x,Y) . 
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III. Clebsch (Legendre) Condition . At each element (t ; x>^jL) of 
a minimizing arc C the inequality 

m 


must hold for every set (Y lV . . . ,Y ) . 

IV. Second Order Condition . The second variation of h along a 
minimizing arc C is non-negative for every variation of C satisfying 
the equations of variation. 

(Cf. Ref. 6, pp l6.) Wo use of this condition is made in this paper. 


YZ H Y.Y. < 0 

TA y i y j 1 J = 
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DEFINITION OF SYMBOLS 


R 

r 

V 
v 

L 

f.g.f.g 

t 

V 

k 

T 

m 

X, y, a 


Vehicle position vector 
Distance to vehicle 
Velocity vector of vehicle 
Speed of vehicle 

Perturbation displacement vector 
Coordinate functions 
Mass parameter 
Time 

Time at which the natural end point is reached 
Magnitude of thrust 
Direction of thrust 
Mass of vehicle 

Lagrange multipliers or adjoint variables 


a Semi major axis 

n Mean motion 


d i 

0 Incremental eccentric anomaly 

L,f 0 ,f„,f. Functions of 0 defined by Eqs. (48) 

(X } Adjoint variables defined by Eq. (18) 

{r } State variables defined by Eq. (18) 

{6(t_) } Residual vector defined by Eq. (19) 

{a} Variational parameters 

{p }, {q } Defined by Eqs. (21), (22), (23) 

[$] Partial derivatives of state variables as defined by Eq. (25) 
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[ A ] Partial derivatives of adjoint variables as defined by Eq. (26) 

[F],[G],[J] Defined by Eqs. (27), (28), (29) 

Subscripts 

u 
o 
E 

A, B 

Superscripts 

k Value at the kth iteration 


Unperturbed solution 
Value at the initial time t 

o 

Value at the natural end point 

Values corresponding to variational parameter set A or B 
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By Jack Richman 


SUMMARY 

/t ?<> 3 

This report contains a description for the solution of the two-point 
boundary- value problem of the calculus of variations for optimum orbits. 

The method employed uses Lagrange multipliers and Pontryagin's 
maximum principle to obtain the decision functions. 

In addition, two differential correction schemes are described. The 
first scheme is a "method by forward integration, " and the second is an 
alternate "method by backward integration" that attempts to reduce the 
difficulties that might be encountered in inverting a differential correction 
matrix. 

The optimum orbit is determined by a perturbation method similar to 
that of Encke and accommodates hyperbolic as well as elliptic orbits. The 
equations necessary for the generation of a digital-computer program are 
derived. 
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INTRODUCTION 


The usual methods of solving the two-point boundary-value problem of 
the calculus of variations involve the use of iterative gradient techniques. With 
these methods, the desired solution is reached only after making a great 
number of incremental variations and examining the changes that these varia- 
tions cause. As one might expect, the rate of convergence for this method is 
very slow. 

Another method of solving the two point boundary value problem of the 
calculus of variations, which will be described in this report, is one where all 
the decision functions and trajectories that are being used are extremals. This 
method uses, in addition to the state variables Lagrange multipliers or 
adjoint variables that play the key role in deciding the optimal direction of 
thrust, time of thrust duration, etc. The adjoint variables also define the 
natural end-point conditions by which the two-point boundary-value problem 
can be terminated. This natural end point, in general, will not be the desired 
end point. A differential correction scheme provides the means of obtaining 
another optimum trajectory the natural end point of which will be closer to the 
desired end point. 
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EQUATIONS OF MOTION 


In a Newtonian system, the equations of motion of a particle that is in the 
gravitational field of N attracting bodies and is subject to other accelerations, 
such as thrust, drag, oblateness, radiation pressure, etc., are given by 




K=1 



— VB 

K 


3 

VB 



(1) 


The problem that will be considered here is one in which the vehicle is in the 
gravitational field of only one body and is subjected to a variable thrust k_. In 
this case, Eq. (1) is reduced to 


R=-M-^+^-T ( 2 ) 

3 m — 
r 

where T is a unit vector in the direction of thrust. The magnitude of the thrust 
is taken" to be proportional to the mass flow and is given by 


k - - cm ( 2 ) 

The constant of proportionality c is related to the more commonly used constant 

specific impulse I by 
sp 


DERIVATION OF OPTIMIZATION EQUATIONS 


In the derivation of the optimization equations, it will be assumed that the 
vehicle can have two possible values of thrust, either k = k max or k = k min . The 

magnitudes of these two thrust values may differ with each stage. 


Minimum -Fuel Condition 


The value of the integral to be minimized is given by 





- ihdt 


I 


( 5 ) 


and the conditions of constraint are given by 


UR k 

V + ~~ - T = 0 
” 3 m — 

r 


R - V = 0 


• k A 
m + — = 0 
c 


Because these conditions of constraint are satisfied at every point on the trajectory, 
we may rewrite Eq. (5), without changing its value, as 


T • jU R k • u 1 

! = J |-*+X -(V+-f- -- T) +y - <R-V)+a (m + f) j dt 

t r 

o 


F L (R, R, V, V, m, m, X, y, a) dt 


where X(t), ^(t), and a(t) are undetermined Lagrange multipliers that are chosen 
so as to determine the optimum decision functions required to solve the problem. 

Applying the Euler Lagrange equation 

A /A, 

dt ^ 5 q . ^ o q . ^ ^ 

l 

to the state Variables, results in the following set of equations: 

X +% ~ 0 


jU X 3/i (R • X ) 

■y - t 7 r-=— R - 0 

" %J O 

r r 


X • T = 0 



:J7 


Equations (6) and (9) can be combined to form 
lill k 

R= =- + — T 

— 3 m — 


u A 3u (R • X) 

i = 'T ■’ T~— ' B 


k - - c m 
k 


a 


2 I' I 


111 


In addition, the natural boundary conditions are 

it„ 


X • 6V 


= 0 


(10) 


y_ • 6 R 



(a -l)6m 



(11) 


Because variations in the position and velocity at the end points are zero, 
the first two expressions of Eq. (11) yield no additional information about the 
values of X and y_ at the end points. The variation of mass at the final end 
point, however, is not zero, i.e. , 6m (tp) ^ 0 . Hence, the only way to satisfy 
the third expression of Eq. (11) is to demand that 

a (t p ) -1 = 0 (12) 

The only additional information that is necessary to completely define the extremal 
is the determination of the optimum thrust vector and the duration of this thrust. 

For the determination of this decision function, we make use of Pontryagin's 
"Maximum Principle, " ("'"’^which states that a necessary condition for an integral 
of the form of Eq. (7) to be minimized is that the Hamiltonian be a maximum. The 
Hamiltonian for this problem is given by 
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H = X • V +Z ■ V + cnh 




R + 


- A * R - 

m — j - - 


crk 


( 13 ) 


r ur-X • • i 


+ k! 


r X • T a 1 


m 


For H to be a maximum, the unit thrust vector T must be in the direction 
of X , or 

X 

T = (14) 

1*1 

Therefore, the coefficient of k in Eq. (13), which is defined as the switch 
function, becomes 


The necessary conditions that must be placed on the magnitude of the thrust for 
H to be a maximum are the following: 


if S > 0 then k = k 

max 

if S < 0 then k = k . 

min 


(15a) 


Furthermore, when thrust is applied, it is desirable to make the switch function 
as large as possible. This can be accomplished by allowing the mass to be as 
small as permissible, which implies the obvious condition that any empty tanks 
or other unnecessary weight be dropped as soon as possible. 

Minimum-Time Condition 


In this case, the value of the integral to be minimized is given by 


I = 


r t F 


dt 


t 

o 


Xpr MR k ‘ k ~ 

J I 1 +x -(Y + -f- --T) +V -(R-V) +CT(m +-)jdt 


(16) 


Application of the Euler Lagrange equations and Pontryagin's Principle lead to 
the exact same results as the minimum-fuel condition, with the exception of one 
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of the natural-boundary conditions. In place of the third expression in Eq. (11), 
we now have 


or 


cr6m 



It 


o 


= 0 


o (t F > = 0 


(17) 


Therefore, for the "minimum -time" condition the natural end point occurs 
when cr = 0. 


ITERATION SCHEME 


General Procedure 


The problem is to generate a set of initial adjoint variables such that an 
optimum orbit can be computed where the natural end point matches the desired 
end point. (The end points are, of course, given by terminal values of the state 
variables.) With initial values of the state variables specified and an estimate 
for the initial values of the adjoint variables, an iterative method can be used to 
solve this problem. Improved estimates for the initial values of the adjoint 
variables can be obtained by computing the residuals or differences between the 
values of the state variables at the desired end point and the natural end point 
and then applying a differential correction matrix to these residuals. We define 
the {r} , 


{X }, and {6 (t„) } vectors as 


{r} = ( 


m 


and 


{fi(t F )}= 




x(t F ) -x E 

y(t F ) - ye 

z (tp ) ~ z p 

* <V - *E 
y (tp) - $e 

z (tp) - Z E 
m (tp) - m E 


.X 

-*y 

-x„ 


U} = { 




J 


(18) 


(19) 


J 
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where the subscript E denotes the values of the state variables at the desired 
end point. 

The Kth approximation to {X (t 0 ) } is designated by {A^(t )}, and it is 

desired to obtain an improved value of {A^ + ^(tp)} . The procedure is as 

follows: using {A'*^(t o )} in the integration scheme, the position, velocity 

and mass at time t^,, as well as the residuals [6^(1^)}, are computed; and 

the initial values of the adjoint variables are then changed so as to reduce the 
residuals. 


(2o> 

where {AA^(t Q ) } is to be found by using a differential correction matrix. 

Methods for Obtaining the Differential Correction Matrix 

Making use of Eqs. (14) and (18), the first two expressions of Eq. (10) can 
be written as follows: 


{i* } = {q ( {r}, {A }) j- or = q ± (C r } , {A } ) 
[A }= {p ( [r], {A})} or A. = p. ({r} , {A }) 


( 21 ) 


where 


M 1 

^2 


= x 

= y 

-- z 


A 


<L. = " 


d* = 


UX 

3 

r 

i 

3 iw* 

X 

1 A 1 

(22) 

_jdX 

3 

r 

k 

+ 

m 

A 

_X 

|a | 


U z 

k 

A 

Z 


3 

r 

+ m 

|a I 



- — 5- + ^r 


k 

c 
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and 



3jU (R • X) 

5 

r 



3/i (R • X ) 



y 



3u (R • X) 

5 

r 


P 4 



(23) 



P 


6 




z 


P 7 


k 


2 

m 


IM 


Taking the variations of Eq. (21) with respect to a set of parameters 
{a} = (a v a 2 , a ? ), we find that 


£ [#] =[F ] [<f] + [ G ] [A] 


^ [ A ] = -[ F ]* [A] + [J ] [4 ] 


( 24 ) 
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where 


d x(t) dx(t) gx(t) dx(t) dx(t) ^x(t) dx(t) 


8a, 


da. 


da. 


da. 


9a 5 ^6 


da r 


^y(t) 

da^ 


dz(t) 

da. 


<•>!$ 


dx(t) 

da. 


( 25 ) 


s y(t) 

da x 


d z(t) 
da. 


d m(t) 
da^ 


d m(t) 
da„ 
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Method of Forward Integration . Two convenient sets of parameters to 
work with are the sets that consist of the initial values of the state variables 
and adjoint variables, which are, respectively 

<“>A = { r(t o>} 

and (29) 

£“ 3 b = {>• <t 0 >} 

Using these sets of parameters, Eq. (24) can be integrated "forward" 
simultaneously with equations of motion, using the initial values of [$>] and 
[ A ] as given by 



The differential corrections are obtained by solving the system of equations 

{ Ar <V} =r* A <‘F>] { Ar <‘o>MVf>] 

( 31 ) 

{AX(t F )} = [A A (t F )] {arty} + [A B <t p )] {AXiy} 

and, because 

{ Ar (t Q )} = 0 and { Ar (t F )} = { 6 (t p ) } 

we find 

Wo>HVf>] { 6 < l F>} < 32 > 

An interesting feature of this differential correction scheme is a tendency 

for the inverse of the differential correction matrix [$ (t_)] to become more 

B F 

and more singular as the time arc increases. This tendency toward singularity 
is a problem of utmost interest. 

Method of Backward Integration . If the use of double-precision tech- 
niques fails to provide the required numerical accuracy for the inverse of the 
matrix, an alternate method of generating the differential correction matrix 
can be used. This alternate scheme employs a method of "backward" inte- 
gration to provide a differential correction matrix consisting of the sum of two 
matrices, only one of which requires inversion to produce the differential 
corrections, In this case, the two sets of parameters consist of the final 
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values of the state variables and adjoint variables, which are, re- 
spectively 

{<*) A = { r <y} 

and , > 

(a 3 B = { x ftp) } 


Using these sets of parameters, the variational Eq, (24) can be inte- 
grated "backward. " The procedure is as follows: the equations of motion are 
integrated "forward” until the natural end point is reached; the residuals are 
computed; and, then, Eq. (24), together with the equations of motion, are 
simultaneously integrated "backward" starting at time t^, and ending at time 


t , using for initial values of [$] and [A]: 



a 


r 

>A<*F> ! 

= I 


L*b <V ] - 0 



S and / 

V 

r / 

lA <‘f>>° 

> 


A <V J = I J 


(34) 


The differential corrections are obtained by solving the equations 


{ 

{ 


i Ar <y} + [vv] { AX( v 

Ar(t o ) } = [W] { Ar ( V" tV'oO { AX(t F>} 


(35) 


and, because, in this case, 

{ A r (t Q ) } = 0 and {a r (Ip) } = {d tt p )} 


solving Eq. (35) for {AX(t Q ) } , we find that 

{£X(t o) } =[CA A (t 0 )]-[A B (t 0 )][* B <to)]' 1 [' I 'A( t 0 ) ] ]{ 6 <tF>} < 36 > 


Convergence of Iteration 

Several difficulties are connected with the above iteration scheme, and 
some of them might be crucial enough to cause divergence of the iteration. 
These difficulties might arise for the following reasons: 

1. In the variational equations, the variation of burning time 
is not accounted for. 
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2. The inversion of a matrix is required in both methods to 

obtain the differential-correction matrix. Furthermore 
this inversion becomes more involved since the residual 
m (bp)~ m g of the vector is unspecified and requires 

additional computation. 

3. The change Atp in the final time has not been taken into 
account. However, this should be included by considering 
the additional transversality condition which results in 
{A.} • {r} = 0 . 


DIGITAL PROGRAM 


Trajectory Equations 

The equations that completely define the trajectory have been described 
previously. The order in which these equations are programmed for the 
general case (with thrust) is as follows: 


I * 


m 


b 


>0 
< 0 


k = k 

max 
k = k . 


mm 


m 


a = 


. is 

c 

is 1 a 1 
2 
m 


^ [<*]= [F][$] + [G] [A] 

^ [A] - - [F ]*[ A] + [Jim 


R = 


jUR k 

r 3 + m 1^1 


X 



3/i(A-R) R 


4 +At 

m =m(t)+ mdt 


a =a(t)-t 


„t + At 


ff dt 


(37) 
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These equations are integrated until the natural end point is reached. At that 
time, the residuals are computed and compared to a predetermined set of 
maximum permissible values {f } . 

If 6 ft ) £ e for all the residuals, then that trajectory is the solution 
j F 7 j 

to the two-point boundary -value problem. If 6 j( t F ) > ^ for any of the re- 
siduals, then a differential correction is applied to the initial values of the 
adjoint variables as described previously. If the alternate differential correction 
scheme is used, then a "backward" integration is necessary before any correc- 
tions can be applied. 

Numerical Procedures 

The differential equations of Eq. (37) can be integrated numerically with 
a Runge-Kutta fourth-order method. To reduce any accumulation of error 
that might result from a number of step-by-step integration, however, it is 
convenient to write the equation of motion for the high thrust case in the form 

R = R u + £ (38a) 


The velocity and position vectors can be written as 




R = R u + L 


(38b) 


where R u is the unperturbed solution and £ is the perturbation. 
In this method, R u is taken as 


R 

= — T 



cm 

(39) 

— u 

m — i 


m — i 

i = 

MR + 
— 

k_ 
m . 

T - T. 

(40) 


r 


Eq. (40) is integrated numerically, and the solution to Eq. (39) is 


R u = f R(t.) + gR (t.) + hT(t.) 
R u = f R(t.) + feR(t.) +hT(t.) 


( 41 ) 
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where 


f = 1 


g = t - t. 


h - “ c {^ [ m lo S m " m. log m. - (m-m.jJ-tt-t.) log m.j 


f = 0 
g = 1 

h = - c(log m - log m.) 


m = m. + (t - t.) rti 


(the subscript i refers to values at time t.) 


This perturbation method, or Encke scheme as it is commonly called, 
will reduce inaccuracies occurring in numerical integration, provided that the 
perturbation terms are small compared with the total solution. Whenever 
these perturbations become too large, a rectification takes place, i.e., an 
initialization occurs in which the values of the variable at time t now becomes 
the values of the variable at time t.. A rectification takes place whenever any 
of the following conditions occur: 


L 

r 


> € 

pos 


(position rectification) 


i 


> <r 


vel 


(velocity rectification) 


( 42 ) 


,[ 


2 T 


T. > e 
— i acc 


(acceleration rectification) 


SOLUTION OF EQUATIONS FOR THE COASTING STAGES 


The solution of the equations of motion and the Euler Lagrange equations 
can be derived in closed form for the coasting period. In the no thrust region 
(k=0), the equation of motion reduces to 


R = - 




3 

r 


(Kepler problem) 


( 43 ) 



The two-body orbit that results from the solution of Eq. (43) with the initial 
conditions 


51 


R(t.) = R. 

R (t.) = R. 

can be written as a linear combination of R. and IT as 


(44) 


R-f R. + gRi 
R=f R. +gR. 


(45) 


The coefficients f, g, f , and g are obtained as follows: we represent the initial 
conditions by the set of elements 



1/2 

u 

= a 3/2 


(elliptic) 


(46) 


n = 


.1/2 


(" a ) 


372 


(hyperbolic) 


This results in the following Kepler's equation 

r. d. 

n(t - t.) = 0 - sin 0 + — sin 0 + — (1 - cos 9) (elliptic) 

1 a M 

r. d. 

n(t - t.) = sinh 0 - 0 - — sinh 0 +— — (cosh 9-1) (hyperbolic) 
x l' a } 


( 47 ) 


where 0(t) is the incremental eccentric anomaly E-E.; the functions f^, f 2 
f , f are defined as 

O /r 
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^(0) = 0 - sin 0 

fg (9) = 1 - COS 6 

f 3 = sin 0 = 0 - ^(0) 
f 4 = cos 0 = 1 - f 2 (0) 

f^(0) = sinh 0-0 
f 9 (0) = cosh 0-1 
f 3 (0)= sinh 9=0 + f ^ (0 ) 
f.(0) = cosh 0=1+ f (0) 


(elliptic) 


(hyperbolic) 


(48) 


and the solution of the two-body problem for both elliptic and hyperbolic orbits 
is given by 

| a | 

f = f 0 + 1 

r. 2 
i 


* - -T f i + - b 


al VlStT ■ 


d. 

i 


v 


/ u a 



1 k 

/ — f | 

1 

f = - 



a 

r. 


* i 

i 

g = - 

M 

f 2 + 1 

r 


r. 

i 


d. 


n (t ~ t i ) - f l + ]aT f 3 


J M & 


(49) 


For the non-thrust case, we also can solve for {A } in closed form. The 
following is a derivation leading to this closed-form solution: the differential 
equation for the adjoint variables are written as 

{A} = - [F]*{A] 


(50) 


where [Fj is defined by Eq. (27); the variational equation for [$] reduces to 


£ [*] = [F][4] 


(51) 


taking the transpose of Eq. (50) and postmultiplying by [ 4] , yields 


ap (xf [«] = - {x}* [Fits: 


(52) 


premultiplying Eq. (51) by {X} , yields 


{X}*^[*]={X3*[F] W 


(53) 


comparing Eqs. (52) and (53), we see that 




or 


dt L 


{X} [4] = 0 


(54) 


Tf. 

Eq. (54) states that {X } [4-] is a constant and, therefore, can be written as 


{ x (t) r [*(tj ] = {x (t K )}* (t K ) j 


(55) 


where t-, is any fixed time in the no-thrust interval; solving Eq. (55) for { X (t) } , 
K 

results in 


{X(t)3 = 



L** (t K ) J { X <‘ K )} 


(56) 


In the case where the set of parameters { a } corresponds to a set of the state 
variables { r] , the matrix [4] can be written as 


0a <*>]= rV*-‘K>] [*A <*K> 


(57) 


taking the transpose and then the inverse of Eq. (57), leads to 
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and combining Eqs. (56) and (58), results in 


U(t)}=[<i>* (t - t K )] _1 {a ( t K ) } 

which is the closed form solution of { A (t) } . 


( 59 ) 


The elements of the j$^(t - t ^)J matrix are obtained by differentiating 
the Kepler orbit elements with respect to R(t^) and R(t^). The elements of 

k “ - i k>- 


5 r «(t) 


A K'Jpq 5r q (t K ) 


, with p, q=l 7, are as follows: 


5x i (t) _ 9x i fJ6 3 |a | f _ 

dx.ft..) Sx . ij + 3 X oj L x i X oi ^ oiJ 

V K' oj J r Q J 


a lx . 

+ — T*-*r* 0 ) 

r 

o 


■3 (t - t R ) +g+— y- (1 -jrr) k 
a 




I p / i i & • C* t 1 1 

— L f„(x. -x ,)x . + f 0 ( + 1 r) — TP x . 

H 2 V i or oj 2 v r Q | a |' 3 Ol 


X 


oj 


(60) 


_ 3x i , 3 la 


dic.(t^) Mi . 
y K' oj 


* V * i ^ojK^oi-^-v^oi. 


| a j ic r 

+ ~ir 1(i ‘i- i oi>r 3 ( t - t K> +g+ pir f 3l 


-x . 


U OJ 


(x.-x .) f 0 +— — f„ X .X . 
' i Ol 7 2 ^r o 2 Ol OJ 



r • L -i 


dx.(t) 3x. . ^|a|x r 

1 “ -=f6.. ;r— -P- I X. +r(x. - X .)— 77 ! —3 (t — t ) + g + -r^T — (1 - 

. ij r 3 r 3 L l ' l oi' |1 J- ' K' 6 I a | n ' |a| ’ 3. 


r r 

o o 


aXj(t K ) sx o . 


|a \k . 
+ 1 1 °J . 


r- I n 


x. + r (x. - x .) — f 0 + -7- (x. - k) x .f 

L 1 '1 or J 2 /i 1 01 °j 


"H f < T7 + lij) X oi X oj + ( *i “ X oi> X oj [ £ + r (1 " ^ f 4] 

o 

(60 continued) 
r • r* r 

[x. + r (x 1 -* oi )^7T^] [ _3(t ‘ *K> + g + T5TC f 3 J 


r o 

o 


dx.(t) M*oj 

= 5x oj ‘ g ij " r 3 


M 


3 f 2 x oj [*i + r (k 1 -*oi)^r-] + /^i-* 0 i>* 0 j[-^ fe+ ¥- f 4- 


oj L r ° r 4. 
J o 


(k. - x )x.f - — fx.x. 
' 1 01' 01 )L 4 01 oj 


where i, j=l, 2 , 3 correspond to the x, y and z components and 

x 0 ' X <‘ K > 

r 0 a r(t K> 
r « r(t) 

The inverse, [4*. (t - t„) ] , can be obtained from the above expression by replacing 

t — -t r -* r 

o 

e - e x - x 

o 

r -* r x — x 

o o 

This results in 



f - g 
g - -g 
f - -f 

g - f 
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. , i>^ \ Summary 

I b 

This report describes an application of a successive 
approximation scheme for optimization of low- thrust trajec- 
tion involving many revolutions about a central body. The 
equations of motion, written in terms of orbital parameters 
rather than position and velocity components, are analyzed, 
and a convenient thrust formula is derived. This formula, 
together with the variation-of-parameters method of trajec- 
tory computation, has been programmed for the IBM 7090. 


INTRODUCTION 

One of the principal difficulties associated with the 
application of a numerical optimization scheme to very low- 
thrust trajectories is the large computer storage required 
when Cartesian coordinates are used. Many points are 
needed to compute each trajectory, and when several solu- 
tions must be stored, the computer capacity can easily be 
exceeded. In addition, the geometry of very low- thrust 
trajectories suggests that the optimum thrust will oscil- 
late in a regular fashion with only very small changes per 
revolution superimposed on this oscillation. We are con- 
cerned with these small changes, and it seems likely that 
they will become obscured over many orbital periods. 
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A natural solution to the problem of storage is to 
employ a variation-of-parameters integration routine. Fur- 
thermore, when the entire problem is rewritten in terms of 
parameters, a thrust formula suggests itself which is a 
function of the parameters and a set of slowly varying con- 
trol functions. It is anticipated that these control func- 
tions will be sensitive to the small cyclic changes in 
thrust . 

The specific problem to which we have directed our 
attention is that of minimizing transfer time for a two- 
dimensional very low- thrust trajectory when only the thrust 
direction is variable. The formula for the resulting 
thrust direction contains three control functions, and 
specific examples indicate that these have the desirable 
properties we are seeking. The following contains the de- 
velopment of the thrust formula and specific equations for 
the parameters we have chosen. 


DISCUSSION 

Let p£ be the parameters which describe the trajec- 
tory at any time t, and let 


dp . 

= ULG i (p 1 , . . ., p n , t, ip) 


i = 1, . . . , n (1) 


be their differential equations where p is a small quan- 
tity, in our case the thrust acceleration, and ip is the 
angle of the thrust vector. If we write the Euler-Lagrange 
equations for minimizing time with ip as the control func- 
tion and these differential equations as side conditions, 
the equations for the Lagrange multipliers, become 



d 


n 


I 

k=l 


9 £ 

k dp. 


i = 1 


y • • • y 


n 


( 2 ) 
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0 


dG, 

•* I \ w ■ 

k=l 


( 3 ) 


Because for our problem if/ appears only in linear combina- 
tions of sines and cosines in the G^, the last equation 
yields a formula for the tangent of if/ in terms of the 
and the parameters p-^. But the p^ are slowly varying, 
and if the partials dGj^/dp^ are not too large, the £j_ 
vary slowly as well, since their time derivatives also con- 
tain the small quantity [i as a factor. Therefore, if we 
adopt the as new control functions using the tangent 

formula to replace if/, we will have substituted n slowly 
varying quantities for one rapidly changing one. Under 
favorable circumstances, this would mean both increased 
accuracy and reduced storage space requirements. 

In practice, it is not necessary that all the parame- 
ters be constants of the zero thrust or unperturbed motion, 
but only that they be slowly varying in some sense and that 
dGk/dpi be of moderate size. In applying our analysis to 
very low- thrust trajectories, we have chosen, for sake of 
computational ease, to integrate directly for 0, the 
angle of the vehicle in the plane of motion, rather than to 
calculate a time parameter such as the time of perigee. An 
examination of the differential equation for 0 shows that 
it is well-behaved and does not contribute significantly 
more to the time rate of change of the than the other 

parameters. In terms of the classical orbital elements our 
parameters are 

p-^ = J a (1 - e 2 ) /k* = h 


p 2 = e cos oj = q 
p^ = e sin oj = s 


( 4 ) 


P4 = d 
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where k = Gm, the universal gravitational constant multi 
plied by the mass of the central body. The differential 
equations for these parameters in a two-dimensional very 
low- thrust trajectory are 


= l[h 2 /(l + q cos 0 + s sin G) Jcos ip 


dq 

dt 


[cos 0 + (q + cos 0)/(l+q cos 0+s sin 0 ) ]h cos ip 


+ -[sin 0 ]h sin ip 
m 

(5) 


ds 

dt 


-[sin 6+ (s+ sin 0)/(l+ q cos 0+ s sin 0) ]h cos ip 


- —[cos 0 jh sin ip 


03 

— = [ (l + q cos 0 + s sin 0 ) / (h k) ] 
dt 


where T/m is the thrust acceleration, a small quantity. 
Substituting Eqs . (5) into Eq. (3) yields a thrust direc- 
tion formula, and it should be noted that because f does 
not appear in d0/dt the multiplier ^4 will not appear 
in this formula. Defining two auxiliary quantities, A 
and B, by 

A = HA h/(l + q cos 0 + s sin 0) ] 

+ i [cos 0 + (q + cos 0 ) / (1 + q cos 0 + s sin 0 ) ] 


+ U^fsin 0 + (s + sin 0)/(l + q cos 0 + s sin 0) ] 
B = H 2 sin 0-^3 cos 0 , 
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we have for tan ip 


tan ip = A/B 


With proper attention to signs, we can solve this equation 
for the sine and cosine of ^ and substitute these into 
Eqs. (5) . When this is done, we have a set of differen- 
tial equations of the following form: 


dp i 

dt~ = H i^ p i' * * * ' p 4* £ l f 



with the as new control functions. 

A successive approximation scheme employing these 
equations has been programmed and check runs have been 
made. The early runs indicate that the estimates of the 
time rate of change of the are valid and that the pro- 

gram can be applied without exceeding available computer 
storage to the problems for which it was designed. Sub- 
sequent efforts will be aimed at assessing the scheme's 
merits and at the possibility of refining it. 
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SUMMARY 


The three dimensional optimum trajectory relations developed by 
Messrs J. G. Cox and W. A. Shaw in Reference 1, are transformed into 
a form that appears more amenable to low thrust trajectory calcula- 
tions. Orbital element coordinates, commonly used in Celestial 
Mechanics, are employed due to their slow variation in low thrust 
applications. Combinations of these elements and a generalized ec- 
centric anomaly are utilized in arranging the resulting equations 
in a form which does not contain circular singularities. 
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LIST OF SYMBOLS 
Semi-major axis 

Abbreviated notation for cosine 

Generalized eccentric anomaly (see Fig. 3) 

Parameter linking the eccentric anomaly with the 
generalized eccentric anomaly E (see Fig. 3) 

Numerical eccentricity 

Thurst 

Gravitational force 
Gravitational constant 
Specific angular mementum vector 

Components of angular momentum in equatorial axis 
system. 

Inclination angle of the plane of motion 

Specific energy constant 

Semi-latus rectum (orbital parameter) 

Mass Of attracting body 

A vector perpendicular to the line-of-nodes vector N 
Vehicle mass 

A vector along the line of nodes directed toward the 
ascending node 

Mean motion 

2 2 2 L 

(q. + q^ + ) 2 y Absolute value of Lagrange 

multiplier vector q 

Lagrange multiplier vector in equatorial coordinate 
system 
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Lagrange multiplier components in equatorial coordi- 
nate system 

Perturbation attraction force 


qp q 2 » q 3 


R 

r 


S 


, 2 , 2 
(x , + X + X 

pl_ p2 p3 
tor x~ 




Absolute value of position vec- 


Abbreviated notation for sine 


t 


Time 


u 

Up u 2 , u 3 

V 


y 


X 

p 


Position vector in equatorial coordinate system 
Equatorial coordinates 

Position vector in the plane of motion, x, y, coordi- 
nate system 

Coordinates of an axis system in the plane of motion 
with the x - axis directed toward the ascending node 
and the y - axis 90 degrees forward 

Position vector in plumb-line coordinate system 


x 


pi’ p2> p3 


Plumb-line coordinates 


Greek Symbols 


«c 


0 

e 

* 

e 


C 

6 

A 


Ai, A 2j A 3 


Notation parameter representing either u or q 
Notation parameter representing either 
Mean longitude at epoch 
e - ft , ,r Mean argument" at epoch 
(1-e^)^ cos i 

Position angle from ascending node 

Lagrange multiplier vector in plumb-line coordinate 
system 

Lagrange multiplier components in plumb-line coordi- 
nate system 
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v 


£ 

a 

n 

0 ) 


True anomaly from perigee 
Eccentric anomaly 
e sin 0) 
e cos u) 

Longitude of the ascending node 
Argument of peri-apsis 


Subscripts 

N 


Parameter is evaluated at the time of passage through 
the ascending node 

Indicates initial condition 

Indicates physical parameter (as apposed to Lagrange 
multiplier parameter) 


A 


Lagrange multiplier parameter 

1,2 


The subscripted parameter pertains 
respectively. 

1,2,3 


Orthogonal cartesian coordinates 

Other Notations 

(’> 


Denotes the first time derivative 

f) 


Denotes the second time derivative 

( ) X 

( ) 

Multiplication 

(") 


Denotes a vector 

(") X 

(') 

Vector product 

(') • 

(") 

Scalar product 


Absolute value 
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INTRODUCTION 


The ’’minimum fuel’ 1 trajectory equations were derived in Reference 
1 by Messrs. J. G. Cox and W. A. Shaw utilizing the following assump- 
tions : 

1. Spherical earth. 

2. Inverse gravity law, F = - 

& 2 
r 

3. The only forces acting on the vehicle are thrust and gravity. 

4. Rotation effects on the rocket are ignored, 

3. Constant fuel burning rate. 

6. The center of mass of the vehicle is fixed with respect to 
the vehicle. 


The equations derived under these conditions are: 


x 


P 




X 


X 



3GM 




x„ 


(i) 


where x is the position vector in the plumbline coordinate system 
described in Reference 2 and A is the corresponding Lagrange multi- 
plier. The computational method of Reference 1 is that_of approximate 
integration starting with initial conditions on x^, x ,A and 


The computational scheme of Reference 1 does not readily lend it- 
self to rapid trajectory calculations for situations in which the 
thrust force is small in comparison with the vehicle mass. This is 
due to the large coordinate variations required for relatively small 
displacements of the vehicle away from the attracting body. Conse- 
quently, expressing the equations in a coordinate system whose natural 
properties are better suited to approximate integration would be ad- 
vantageous. This study is then a exploratory investigation into the 
possibilities afforded by the elliptical element coordinate systems 
used in Celestial Mechanics. 



69 


COMPARISON WITH THREE BODY EQUATIONS 

The three-body perturbation equations (Reference 3 or 4) may be 
obtained in the form 


o , , x, 

T + k (M + V -7=13 = 

X 1 [xf 

o X 2 

“ + k Z (M + m 2 ) — p = 

x 2 Uo 


dk 1.2 
d x. 


where and are the position vectors of the bodies of mass and 
1^2 respectively and the * (i = 1,2; j = 1,2) are the perturbation 
attractions of body i on body j . 


Equations 2 are usually referred to equatorial or ecliptic axis 
systems and it is convenient to transform Equations 1 into an equa- 
torial system for the purpose of comparison. The transform relations 
(Reference 2) are: 

Mi K - 90 l2 ^ 


* "[♦<>], k 


- 90 °j A 


where u is the position vector and q is the Lagrange multiplier in 
the equatorial system. The "minimum fuel" trajectory equations in 
the equatorial system then become 


L 

+ 

GM- 

u 

= F 

u 

k I 

mp 

A& 

+ 

f GM ' 

q 

_ 3GM 

q 


3 

1 r J 


r5 


(u . q)u 


where p = |q | and r = |u | . Adding =EL. q to both sides of the second 
of Equations 4 yields 


u + r 3 u mp q 


^ T 
P , 


q = - GM 


h> 9 - ^(u-q)i 

rJ 
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A comparison of Equations 5 with Equations 2, neglecting nu and 
with respect to M, indicates they are of the same form with = u, 
x 2 = q, k 2 M = GM,3 2^ x 1 = <* R r / 9 ^ = 3 F/mp, and 3R 2 ^3x 2 = 3 R^/3q = 

- GM £(l/r 3 - 1/p 3 ) k - 3/r 5 (u • q) u] . 

The classical variation of parameter procedure, described in Refer- 
ences 3 and 4, may thus be adapted to the trajectory equations, Equa- 
tions 5. 


ORBITAL ELEMENT EQUATIONS 

The application of the variation of parameter method to Equa- 
tions (2) is given in detail in References 3 and 4. The procedure 
results in two sets of six first order equations in six orbital para- 
meters. The form chosen here has a - the semi-major axis, e - the 
eccentricity, e- the mean longitude at epoch, i - the inclination 
of the orbit plane, cu- the argument of perihelion and ft- the longi- 
tude of the ascending node as variables. Only one set of equations 
is presented; however, it must be remembered that this set actually 
represents two sets of the same form, one representing the real or 
physical situation coming from the first of Equations 5, and the 
second representing the Lagrange multiplier equations coming from the 
second of Equations 5. 

The equations are 

a = -2- -as 
na 3 e 



a 
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Equations 6 contain singularities at e = o and e = 1, (i.e. for 

both circular and parabolic orbits), and as such present computational 
difficulties for near circular or near parabolic trajectories. Sing- 
ularities and the corresponding computational problems are also pre- 
sent when the plane of motion is in or near the equatorial plane. 
Consequently, Equations 6 are not suitable for trajectory calcula- 
tions. The use of a new set of variables, formed by combining the 
variables used in Equations 6, will allow the expression of the set 
of equations in a form devoid of the circular singularity which is 
most important in low thrust trajectory considerations. 


The new variables which allow the expression of the equations 


with no circular singularities are a , ft , e 


,5 , and £ where 


o = e cos co, £ = cos i 


£ = e sin to , e ,v = e 


(7) 


and a and ft are the same as defined previously. 

The 3R/ 3 (old elements) in terms of 3R/ 3 (new elements) are deter- 
mined by utilizing the chain rule and evaluating the partial deriva- 
tives of the old elements with respect to the new from the relations 
of Equations 7. 


This operations yields ** 

— I = 3 R/ 3a 

3 a 


3 R l 
3 e 

3 R 
dQ 


= 3 r / 3 ( 


3R _ J)R 
312 3e * 


( 8 ) 


3 R 
3 e 


a 

e 


M 

3 a 


+ 


3R 

e 3£ 




JR 


Jr 

3co 


- a 


3 R 

35 


* 3 o 


(l-e^)' s sin i -2— i - 

3 C 




The parameter e is used throughout the remainder of this report 
to indicate e = (a ^ 
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where the parentheses around a partial derivative, (|— &) etc, , indi- 
cate the parameters in Equation 6, The time derivatives of the new 
elements in terms of the old are obtained from Equations 7 as: 

a = (a) 


e * = O - (ll) 

0 - (o) 

a = -f (e) + («) K 

e 

5 = (e) + (oj) o 

C = 6 C ° S - (e) - (i)(l-e 2 )^ sin i 

(1-e 2 )^ 


Substitution of Equations 6, 7 and 8 into Equations 9 then yields the 
desired set of non-singular equations. 


na 


9 R / 3 e * 


e * - 



f L±z*hl 

n + (l-e 2 ) 


_/ 

V 

"2 ' • 




H 


2a jr 

9a _) 


a 


i _m 

na- 3 c; 


_i_ . - (1-e 2 )^ 

9 \ 9 t 

na 1 -f (l-e z ) 2 



(l-e 2 ) 


h 


8R 

35 ' 


-I- ■ * + (1-e 2 )^ 

na“ ; 1 4- (l-e 2 ) 2 



5 



( 10 ) 



73 


The R-^ 2 an( 3 R 2 1 t ^ ie t ^ iree body equations, Equations 2, have 
precise definitions! However, in trajectory calculations only the 
partial derivatives 3R r /3u and 3R^/3q are defined. The a R / 3 (element) 
terms in Equations 10 are determined from these definitions by means 
of the chain rule. Hence, it is necessary to define u and q in terms 
of the elements a, e * , Q y <3 £ and;; . 


The orbital elements of a particular orbit may be defined in 
terms of a rectangular cartesian coordinate system in the plane of 
motion. The procedure used here is a slight deviation from the pro- 
cedure given in References 5 and 6 . A rectangular coordinate system 
may be defined with the x - axis toward the ascending node and the 
y - axis 90 degrees forward in the plane of motion. The radius vec- 
tor r may then be expressed in terms of the orbital parameters a, 0 , 
£ and a generalized eccentric anomaly E which is defined (Ref. 6) 


E 




(in 


where the subscript N denotes nodal passage. Restricting E such that 
its value is zero at nodal passage, the position vector V may be ex- 
pressed 

— ; c 

v = ~ rea ) cos E + v N v N (^_) 2 s in E - ea (12) 


where e is a vector of magnitude e directed toward the perihelion. 
Utilizing the relations (Reference 6 ) 


and 




[ ' 

s 

X 

l 

1 n 1 
! 0 \ 


II 

o 

i 

(T> 

II 

vCX 

l 


v 




(l-e 2 )^ 


l+o 


The posit] 






v = 


X 


2 

(1 ~ -+-) cos E - 
1+0 

(l-e 2 ) % l „ in E ‘ 


1+0 

y 

= a 

5 cos E + ( 1 -e 2 )^ 

sin E - g 

0 


0 



(13) 
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The generalized eccentric anomaly E is connected to the eccentric 
anomaly E usually used, through the relation 

E = E + E* 


as shovn in Figure 3* In the above procedure, E* is chosen to make 

% = °- 


A vector in the x, y, axis system is related to a vector in the 
equatorial system through the direction cosine matrix A* , where 


A* 


s n Cfici 
c n -snci 

0 -Si 


-CftSi 

SflSi 

-Ci 


(14) 


the position vector in the equatorial plane is then expressed 

a = ^A*j 6 (15) 

where a represents either u or q and g either v r or v^. 

The 9 R/ 3 (element) terms then become 


3a 



V- 


-tel 

8 aj 



/ ' 

r : 1 1 


_3R 

= 2_R i 

3 La*] ! 


3 Q 

3 5 1 

L 3n 

6 S 


3 R = 

3 R f 

A* 

1 -il-1 



3 e * 

3a ; 


3 e*J 



3R _ 

3 R 

[a*] 

Jl + 

3 [a*] 

3 i 

3a 

3a • 


3a 

3 i 

. 

3o 





r i 

— i 


-M = 

Li) 
3a i 

Fa* 

] _al + 

[i[a*L 

U 

35 

- • 

1 3? 

3 i 

35 


e y 

6 r 


35 


3_R f 3, A* 
3a 3 i 



6 


} 


(16) 
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The partials of 


A* required in Equations 16 are 


] fs 8 -C^Si -cncil 


= | c fl 3 ^ Si SflCi 

o i 


: 0 ~Ci 


and 


I fc Q 
=Ls a 


3 . 

9 Q 


Si 


- SfiCi SftSi 

-cnci cnsi 

o o 


(17a) 


(17b) 


The partials of 3 required in Equations 16 are 


where 


2 ^ 


l+o 


1 + 0 



a l 

1 £_ = 

a 2 

9 a 

“3 

5 SE - 

o - 


(18a) 


+ o 


- 1)SE 


( 1 ~ e2 ) \ CE 
1 + cj 


Xi 


-SITe fl +n ) - ( 0 + e 2 ) SE - E (1-e 2 ) (l-CE)j l 
2 £l + ( a + e 2 )CE - 5 (1-e 2 )^ Se] 


a 2 = K CE + (1-e 2 )^ SE - C - ^(1-e 2 )^ CE - ? seJ- 


X (l Epl+g) - (g+ e 2 )SE -£ (l-e 2 ) % (l-CE ^I 

[2 [I +a- ( 0 + e 2 )CE - Z (1-e 2 A Se] J 


s = o 


-M = 


ae* 1 + 0 _ ( 0+ e 2 )CE - ? (l-e 2 )^SE 


(C 2 -o-l) SE - £ (1-e 2 )'^ CE 

(1+ °) [(1-e 2 )^ CE - 5 SE] 

0 

(18b) 
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M. = a 

3 a 


where 

Y = £ 2 CE + g (l-e 2 )^SE 

1 (l +°) 2 


>1 


(1 +.g. )| a rS E 


( 1 +a ) ( 1 -e 2 )' 2 


(18c) 

I fe 2 ~ o -l)SE-£ (l-e 2 )^C E 
1 + a 


x -C(l +q ) 2 -g 0(l-e 2 )' 2 SE -£ (1 + o -£ 2 )(1- C E) 

(1 4- a) (1-e"")' 2 jl +a - (a 4- e 2 )CE- (l-e 2 )%SE^J 


2 = 


a - SE - + £(l-e 2 P CE- £ SE J, 


(1-e ) 


J 


v ( C(l + o 'll -€ 2 Ja-e 2 )^SE -£(l+o - £ 2 ) (1-CE) \ 
L (1 +a ) (1-e 2 )^ |T + a - ( o+ e 2 )CE- (l-e 2 )^SE ]J 


3 = 0 


where 


M = a 

3? 


«1 

6 2 

-S3 


(1 +a)(i- e z ) 


-hi 

1 + O 


(18d) 


6 L = L SE (l-e 2 ) i SE + 2ECE + C ( ^-a -l)SE-£ fl-e 2 ^f!K 


1 + ° 


-} 


vf.(l-e 2 -g 2 ) (1-CE) +2 £ (1-e 2 )^ SE 

l (l-e 2 )^ [l +0 - (a + e 2 )CE - £ (l-e 2 )^SEp 

60 = (CE-1) - -L-SE_ + f(i_ e 2 )^ CE - £Se1 

( 1 -e 2 )^ J 


X q-e 2 -r 2 ) (1-CE) +2 £ (l-e 2 V 2 SE 


(1-e 2 )^ fl + a- ( a + e 2 )CE - K (l-e^SE, 


r 7 } 

2 seJJ 


0 
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Finally, the partials of i required in Equations 16 are 


3i 

3a 





(19) 



1 

(l-e^)'2 sin i 


The two sets of SR/ 3 (element) terms in Equations 10 are 
then 
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and 



ULi 



9 _% 

9a x 



3R a 

9 ^X 


- GM 




- GM 


l ( T 

r -> 




5> j ] i** ]J»i; s 3 

5) "]L'n7 M]°x 


- GM 








- GM 




_3 

r5 


(H . q) 




The partial derivatives on the right hand side of Equations 20a and 
20b are evaluated from Equations 17, 18 and 19. 

The complete form of the equations of motion could be formulated 
by the substitution of Equations 20a and 20b, evaluated through Equa- 
tions 17, 18 and 19, into the respective sets of Equations 10. 

Space considerations preclude presentation of the complete form. It 
may be noted however, that the only singularities in Equations 10 are 
contained in the 3 R/ 3 (element) terms. An examination of the right 
hand sides of Equations 20, in conjunction with Equations 13, 17, 18 
and 19, indicates that 3R r /3a r > 3 R r^5r anci 3 R r/3£r contain singulari- 
ties at i r = 0 and e r = 1 and that 3 / 3a * 3 R^/3 £ A and 3 R a /3 £ A contain 

singularities at i^ = 0 and e A = 1. However 9 the circular singulari- 
ties, e r = 0 and e\ = 0, present in Equations 6 have been removed 
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through the transformation. This form of the equations should then 
have some usefulness in the calculation of low thrust trajectories 
that traverse orbits which are circular or elliptical and do not lie 
in the equatorial plane. 
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PRELOAD RELATIONS 


Initial conditions for a trajectory calculation will usually be 
given in the plumbline domain. It will therefore be necessary to per- 
form the following preload calculations to specify the initial condi- 
tions in terms of the orbital elements, 

1. Position and velocity vectors in equatorial axis system: 


u o = 


^ °J i [A 0 " ^ Xp 0 


“o - [*oj 1 [*o - 90 ] 2 ^ 


po 


Lagrange multipliers in equatorial axis system: 

% - [*~j l fo - 9 °°] 2 *„ 

% = [♦ 1 [Ao- 90°J 2 t 0 

2, Specific angular momentum, H 

IK = u xu 
r o o o 

Ko = q o x % 

3, Specific energy constant, k 

k ro = T % 2 - 


k 


X o 



GM/q 0 



4, Semi-major axis, a 


a ro " - » l/2k ro 
a Xo * - “ /2k A o 

5. Semi-latus rectum (orbital parameter), £ 

‘ro * H ro 2 /GM 

*Xo * \ 0 2 /™ 


6. Eccentricity, e 


^x*o ” ro / a ro) 


6 Jo - (1 ‘‘jo /a xo^ 


7. Eccentric anomaly, 


-1 

7? — p A C 

^ a ro 

- u o 

3 — tUo 

ro 

a ro 

e ro 

-1 

~ , = cos 

A o 

/ a Ao 

So 

- % 

e \° 


8. True anomaly from perigee, v 


v 


ro 


v Ao 


2 tan 


2 tan 


% 


,1 + e ro N 
( ) tan 

1 " e ro 


1 

T “ 


( 1 + tl2 )% tan i 
1 - e *o 


ro 


A o 
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9. Postion angle from line of nodes, 0 

A vector along the line of nodes in the direction of the ascend- 
ing node, N , is determined by crossing a vector normal to the equa- 
torial reference plane into a vector normal to the plane of motion, H. 


— “ 


— — 



0 


H 1 


H 



1 


2 

0 

X 

H 2 

= 

" H 1 

-1 


h 3 


0 _ 

— — 






A vector perpendicular to the line nodes, M, is determined by cross- 
ing the angular momentum vector, H, into the line of nodes vector N. 


M = H X N = 
o o o 


H 1 H 3 


H 2 H 3 


“ (Hj 2 + H^) 


Taking the scalar product N q with the position vector gives 


'ro * 

u o 

^ro u o 

cos 0 ro = 

u lo 

H 2 ro - 

u 2o 

H lro 

[ Ao • 

^O 

= N Ao q o 

cos 0 Xo = 

q lo 

H 2Ao - 

q 2o 

H 1Ao 


Taking the scalar product of M 0 with the position vector gives 

2 2 

M ro ’ u o = M ro u o sin 0 ro = ( u lo H lro + u 2o H 2ro> H 3ro- u 3o ( H lro +H 2ro) 

q o - ^X ° q o sin 0 X° - ^lAo ^2o h 2X o) h 3Xo -c 13o^ 1X o +H 2xo^ 

from whence 

q __ ^ an -l (ujo ^i ro u 2o ^?ro^3ro ~ u 3o(^-lro ^~2rp) 

( u lo H 2ro - u 2o H lro) H ro 

(qio H 1A o + q 2o H 2 Aq) H 3Aq ' q 3p( H l Xo + H 2 Xq) 

(q lo H 2Xo " q 2o H lXo^x° 


0 Xo = tan_1 



83 


10. Argument of perigee, w 


ro 


= 0 


ro 


ro 


“Ao = 0 Xo ' v Ao 


11. Parameter a 


ro 


= e rQ cos a) ro 


X o 


= e- 


X o 


cos co 


Xo 


12. Parameter £ 


£ ro e ro s '*' n w ro 
“ e X o sin W A° 


13. Inclination angle, i 

Taking the scalar product of the angular momentum vector H with 
a vector normal to the equatorial plane 


H 



H cos i = -H^ 


The scalar product of the vector normal to the line of nodes with a 
vector normal to the equatorial plane gives 


M 


0 

0 

-1 


2 2 

M cos (90-i) = Hj+H 2 


“ro 


= tan 


-1 


i, = tan 
Xo 



From whence 
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14, Parameter, £ 


2 

C _ = (l-e ro ) 2 cos i 


ro 


ro 


n 2 

r = (1-e ) cos i 


* Xo 


Xo 


Xo 


15. Generalized eccentric anomaly, E 


E =2 
ro 


tan -1 a + e2 ro> % tan 7- e ro 

1+a +5 tan -1-6 

ro ro 2 ro 


Ea =2 tan 
A o 


-1 a - e 2 xo) % tan 2 'lla 
1 +a ro + g o tan 2^ e ro 


16. Mean argument from line of nodes at epoch, e* 


e = 0 
ro 


e* = 0 

Xo 


17. Epoch time of nodal passage, t. 


No 

.2 2 . i 

, -j-i , a ro + e ro . ^ 

= t rt - E + sm E 4 


o ro 1+a ro l+o s ro ro 


0~ e rg) r ( 1 - COS E ) 
^ ro ro 


ro 


ro 


a , + ef (1-el ) 2 

= t_ - E. + -^-2 & sin E % + ft— _ (1-cos E Xq ) 


N Xo ° Xq a , Xo Xo 


1 +a 


X o 


18, Longitude of ascending node, ft 


The scalar product of a vector along the u ^ axis with the line 
of nodes yields 


0 

1 

0 


N 


N cos ft = - 
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The scalar product of a vector along the axis with the line of 
nodes yields 

= N cos (90-12) = H 2 




Relations 4, 11, 12, 14, 16 and 18 give the required initial condi- 
tions for calculations in the orbital elements. 
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COMPUTATIONAL METHOD 

Assuming the Runge-Kutta integration procedure will be used, the 
procedure presented below may be used to perform claculations going 
from the jth time step to the (j + 1) time step. The more cumbersome 
relations presented in the previous sections of this report will be 
referred to in this section with the reference followed by an asterisk 

(*)• 


From step j , the following quantities are known: 


a r> 5 ®r ’ 

a r ’ 5 r 5 C r 9 

U l, U ; 

2 ’ U 3’ X P 2 ’ *P 3 * « r’ 


a A> fi A > t.* 

° A> ^ X 5 C A 5 

qj[ > 1 I 3 * ^ ^ ^ 2 * ^ i x * 


i r , e r> [a* 

j r > E r , x r , 

y r > iv 


e A’ 1 A * J X ’ 

V V y x’ 

n 

\ . 
| 

1 


Compute : 

3 [ A * 1 r and 
3 i_ 

si a* 

9 i 

1 ^ (Equation 

17a)* 


r 

r i 

) 

1 ■ 

| 


Compute : 

■jbl-ix and 

a|A* 

lA— (Equation 

17b)* 


3C2 r 

3n x 




3. 

Compute . 

3 v r 
3a 

< 

and - 

) V 

3 a A 

4. 

Compute: 

9 v r 

- and 

3 "a 



3 e* r 

3 e* 

.3. 

Compute : 

3 v r 

o r-j rl _ 

3 "a 



3o r 

cl LIU. 

3° A 

6 . 

Compute : 

3 v r 

and „ 

3v a 



35r 


35 a 


A 





ctg i r 



0 \ 

T ctg 1 A 

1-e 2 A 

A 


(Equation 18a)* 
(Equation 18b)* 

(Equation 18c)* 
(Equation 18d)* 


7. Compute: 
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8. Compute: 


9. Compute: 


3 i r 

? r 

ctg 

i-r 

TcT 

l-e r * 

3 i x _ 

^x 

ctg 

L X 

9 ^X 

w 

3i r 

1 



3 C r 

d-e r 1 2 ) % 

sin 

ir 

3i x 

1 



Hx" ‘ 


sin 

L X 


10. Compute: 


11. Compute: 


v . 


x r 

Yr 

0 


x X 


m-j = mj_! - m (t. 


- Vi) 


12. Compute: 


3R r 3R r SRr 3 R r 3 R,. 9 Rj- 

iZ’ 3 ?*“’ 37 ;’ sr;* 37 ; 


(Equation 20a)* 


13. Compute: jjJVx 9 jfr 9 R X 

3a x ’ 3n x ’ 3£?\’3a x ’ 3? x ’ 3£ x 


(Equation 20b)* 


14. Compute: 


3 R r 


n r a ^ 


3 £ 


'jr 


it. 

e r 


—4- 


f d-e r 2 )" 2 

f 3 R r 

3Rr1 

3 R r 

_i _r . — 

-2a 

1 l+(l-e r 2 )^ 

kw. 

L 0r 3c r ^ r 

3 5rJ 

'Hr 

r 3a r | 


1 9Rr 

n r a r 2 9 
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5r 


1 f ( l-e r 2 )^ o r 9 Rr 


n»-a 


r a r (_ l+(l-e r 9er 

-1 f ( 1 -e r 2 )^ ^ 


7-r- — 5 .- + (l-e r 2 ) * 


9Rr "I 

J 


-1 / (l-e r Z )^ 5 r 9 Rr ^ 3^1 

n r a r 2 ^l+(l-e r 2 )^ 3e* ^ e r ) 3 o r j 


5r = 


_1 f 9R r 

l 9 ^ 


- ? 


n r a r 


3Rr 

r / 


15. Integrate each relation of step 16 to obtain the values of a r , 
Er > ft r > cr r >£ r» and C r at time J + !• 


16. Compute: 


9R i 


^ a A 


9 €* 


f ' / 1 _ 2\ k 


1 f a-o 

n A a^ L W^ 


I£a + a 3_' 

9o x a?x 


9R 

+ C, - 2 a A 

x >t, ^ 


A 9 _^. 


A ' 


*■ A 


-1 


a r, 


n A a A 


3? 


-l foViS 
n x a »t 1+<1 - e » 2) " 

— -i- 


3R > 
8 c x* 
9 R-i 


1 \\ 


+ (l-e A z ) 


n A a A 1+ ( 1 ” e A ) 2 9e * 


9 R. 

Ha j 


) 


(l-e A 2 )^ 


9 R A 
3a . 


■} 


9R 


n A 3 A l 3n A 


£ A ; 


9 R A 
9e * 


} 


17. Integrate each relation of step 16 to obtain values of a^, 

> £ a and c A at time j + !• A 
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18. Compute for time j + 1: 

n r = (GM/a r 3 )% 

n x = (GM/a^ 3 )* 

19. Iterate for E r at time j + 1, from: 

0 r +e r 2 (l-e r 2 )^ 

n r (t-t Nor ) +€* = E sinE r - S 

l+a r 1+ Or 


20 . 


21 . 


Iterate for E x at time j+ 

n A + e \ = E A - 

Compute for time j + 1: 


1, from: 

2 

o,+e/ 

— £ A- sin E x - 

1+0 A 


< 1 - e A 2 )^ 

1+0 X 


e r = ( o r 2 + 5 r 2 )^ 

eA = < 0 A 2 +? A 2 > % 


22. 

Compute 

for time 

j + 1: 


“r = 

_ -1 
tan 

( K r / ° r ) 


1! 

3 

_ -1 
tan 

< 5 a /o a ) 

23. 

Compute 

for time 

j + 1: 



cos -1 

5 r 


L r = 

(l-e r 2 )^ 



COS“l 

^ A 


1 X = 

(l-e x 2 ) >2 


(1-cos Er) 


(1-cos E^ ) 



90 


24. Compute for time j + 1: 

f ? r 2 
Xr ' at [ <u 

X x" a x[ (1 ' 


5 .) cos E r - <Uer 


2 )U r 


sin E, 


l+o r 


£ 2 (l-e^) 2 C^ 

— i ) cos Ex Sin E 

1+0 X 1+0 X 


25. Compute for time j + 1: 

y r = a r cos E r + (l-e r 2 )^ sin E r - 5 r ] 

y\ = a x[ 5 x cos E * + (1 " e x 2)Z sin E x" 5 x] 

26. Compute for time j + 1: 

sjv 


Mr - 


w 


Cfi r 


sn ; 

cu. 


Cllj-Cij- 

-SQ r Ci r 

-Si,- 


Cn X Cl X 

-Sfi x ci x 

"Six 


27. Compute for time j + 1: 


” u l" 


m 

u 2 
L u 3 J 

- Hr 

u 1 

^ o 

1 


28. Compute for time j + 1: 


'll" 

II 

> 

* 

X x 

q 2 

y\ 

Kl 


0 


29. Compute for time j + 1: 


cfl r si r 



" C ft JS i * 
SQ A Si A 
-Ci x 


] 


~ x lp' 



r _ ^ oi 

T 

u i' 

x 2p 

= 

[[ # ° 1 1 

A o - 90 


u 2 

X, 

3p ' 



: 

u 3 
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30. Compute for time j + 1: 



31. Compute for time j + 1: 



33. Compute for time j + 1: 



34. Compute for time j + 1: 



35. Compute for time j + 1: 
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36. Compute for time j + 1: 


fK;. 

! * 3 j 


ejjAo-go®] 


I ?i 

I ?2 

i q 3 


It may be seen from the above procedures, that all quantities 
required for the next time step^, plus the parameters Xp , Xp£> Xp , 

x Pl’ x pn’ X P3 ? 9 ^2’ ^3 ? ^1’ ^ 2 y anc * ^3 are calculated. Z ^ 
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, /L ABSTRACT 

)L^ C 

This paper discusses the problem of the optimum burning program for 
the vertical climb of a rocket. This problem is of engineering interest in 
view of its applicability to the case of sounding rockets or to the study of 
the vertical climb of a rocket prior to pitch-over maneuver. The essential 
objectives of this work, from the standpoint of both physical and variational 
aspects, are: 

I. A study of the optimum burning program based on a generalized 
model. That is, assuming an arbitrary aerodynamic configuration 
and an arbitrary atmospheric scheme. 

II. An analysis of the numerical solution of the boundary- value 
variational problem. In particular, the determination of corner 
points and the integration of an admissible set of adjoint variables 
for the case where the aerodynamic drag is of the form D = D(v, h). 

In this paper a generalized expression for the optimum burning program 
(or control variable program), valid for any arbitrary aerodynamic 
characteristics and any arbitrary atmospheric model assumed, is derived. 

A numerical method for determining the position of the corner points, for 
cases where D = D (v,h), is discussed. Numerical examples are included. 

The case of maximum final altitude is analyzed in the applications 


presented showing the numerical integration of the variable multipliers 
along the extremal, as well as the switching function HJz) . Thus, ; 
practical example of the numerical treatment of the Euler equations and 
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/ b? 6 y 


corner point determination, which is of value as a basic model for the 
understanding of more sophisticated problems not affording closed-form 
solutions, is given. 
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1. equations of motion and variational problem 

The equations of motion, referred to a cartesian system fixed to a 
flat Earth, are written 



( 1 ) 

( 2 ) 

( 3 ) 


Where the following dimensionless variables have been included; 







rn = JH 


The aerodynamic drag of the rocket will be expressed in the following 
general form: 


5 = -kyz*C D (Z,h ) = Sfrj) 


where 


7 - 7 & 


■£ _ / A S 

5 ^ ~ 


( 4 ) 
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It is assumed that the velocity of the gases at the exit section of the nozzle, 

, and the acceleration of gravity, , are constants. Any solution 

of the set of equations (1) to (3), is expressed in terms of the state variables, 
Z(Z) t h fz) J 7n (Z) and the control variable, 

jQ (Z) , The problem under discussion here is that of ^finding the 

optimum solution y ^ ^ ^ 

of Eqs. (1) to (3), satisfying given boundary conditions of the form, 

fa ,n x> 7 I ,i F ~h F ,^ F ,7 F )^0 (5| 

and minimizing a generalized functional of the terminal values 


yj 


G ” Gfe >h x , 7n x , 7 X , , h F , 7n F } ? F ) , 

In the development which follows the generalized state variables will be 
denoted 


( 6 ) 


Cf , i = 3 • Thus, 

T i 


% - * 


% 


= h 


> Q = m 

r 3 


The Euler- Lagrange sum is 

lOj = yi. / > l = 

The canonical variables ( z > i >/> i ) related to 

C 

( T ’% by the equations 


(7) 


t L 


\ - •& r 

i 


( 8 ) 
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are now introduced. Applying the Legendre transformation of the variational 
problem into canonical form 


v c = ^ , nz‘-A-H 

V L 

and forming the Fundamental Function 

d- >, % 


the following canonical equations of the extremals may be derived (Refs. 


1, 5 and 10) 


y. 

L 


) 



= o 


( 11 ) 



( 12 ) 


Eqs. (12) are the equations of motion while Eqs. (11), in this case, are 

identical with the Euler equations since = A* -A , , as derived 

% 

from Eqs. (1) to (3) and (7). Assuming that the control variable is 
bounded, i. e. , ftmin ^ J & — J^max 9 following equations associated 

with different admissible control variations (restricted or one- 

sided admissible control variations and unrestricted or both- sided admissible 


control variations) may also be obtained 




a) " dH 

A 

= o 

A 


A* 7^ " S A«- ’ /- VarMe 



< 0 

for 

S/£0 , 

/ */ 8 »/V ” <ro " 3 ^ 

> (13) 

c) w 
*/ 

> 0 

for 

S/£0, 

~ ^ i 

/3 = X3 - const. 

/ / ^<7*. 

J 

| 

From Eq 

s. (1), 

(2), l 

(3), (7), (9) and (11) it follows that 


X - 

A 

on 

* 

A “ ° 


(14) 

X - 

A 

on 

h 

= o 


(15) 

x + 

A L 

on 

(»■ 

-A) -o 


(16) 

Also, from Eqs 

. (1) to (3) and (9) 




H = 


■ -a)A - A fi * ) *A z 

(17) 


and therefore 
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Eqs. (1) to (3), (13) ( (a), (b) or (c) according to the admissible control 
variations J and (14) to (16) constitute a set of first necessary conditions 
for an extremal. Every admissible solution, Z (z) y h (?) , (T) > 

y3(r) J y U i (T) , (?) > of the preceding se^belongs 

to an extremal. This set is determined since we have 7 equations in 7 
variables. An admissible extremal must satisfy other necessary conditions 
in addition to being a solution of the previous set. We will now consider 
these necessary conditions. 

From the parametric formulation of the variational problem (Refs. 1 
and 5) we have 

-d- f ifii - q’ tfl ] = djQi- (19) 

ch: l c rj dr 

Since in our case ifh = O » it follows that 

r 

(a - + A (4 + j-A* - a - (20) 

along the extremal. Eq. (20) is a consequence of the Euler equations (14) 
to (16) and may replace any one of them if desired. 

1. 1 Weierstrass Condition and Maximality Principle 

The Weierstrass condition requires that at any point on the extremal 

V = Aifb - Ay’ tQ/ , ^ O (21) 

L % 

and since 

Aifl = y Ay’ ~ aH 


(22) 
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then Eq. (21) leads to 


W= J>. Af’ - AH - Af dl. - - AH i O 




Consequently, 


AH = Hfz, <? ,/+a/, ».) - H[r,y ,/*, *) & o 


(23) 


^ is the optimum control. Eq. (23) is the canonical form 


where 

of the Weierstrass condition , and implies that for any admissible set 

^ ^ Q J the optimum control ^ is that which maximizes 

'i 

the Hamiltonian H. This condition, applicable for bounded or unbounded 
control y5 , is known as the Maximality Principle . 

From Eqs. (13) [ (a), (b) and (c) J , (17) and (23) it is seen that in 

the present case the Hamiltonian is linear in the control variable and that 
the conditions previously discussed may be graphically represented as 
indicated in Fig. 1. 

1. 2 Non- Singularity 

An extremal arc of class D* (Ref. 3) will be called non-singular if 
along each sub-arc the determinant 


o2> = 


lDj , 


? ( A 


0 


r 


£ 


V 


, 5 , (24) 
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is different from zero. Then, along each sub-arc, ^ and JbL • will 
be continuously differentiable. For our problem the determinant (24) is 


o 


0 0 0 o 


0 


o 


o o o 


oZ) = 


o 


o 


o 


o 


o 


o 


o o 


o o 


= i (25) 


o 


o 


0 0 4 0 


o 


o 


o 


O O 1 


Thus, any admissible extremal solution obtained will be non- singular . 

1. 3 Transver sality Condition 

In addition to the necessary conditions for an extremal just discussed, 
we can determine from first variation arguments, that at terminal points 
of the extremal, the following sub- conditions of Transver sality must be 


satisfied 
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A G z + A J d *«=° > (*. % 


OC 


(26) 


( d ° S f + A„) dh" = o ,(X e G y a) dt a = o 

In Eqs. (26) the (— ) sign is applied when OtZ = I and the (+) sign 

is applied when 0<r = F . Eq. (26) must be satisfied for any admissible 

S et ( d(f , dy , d 7 I , dZ )? (o , 0 , 0 ,o) 

L 1 C F 

consistent with the equations of terminal variations derived from Eq. (5) as 


(27) 


S ^ d^ y •*£ d^ + ^ dZ j y dZ F = 0,f*Wi7, 

For a normal extremal, as assumed in this case, there exists a unique set 
of variable multipliers. The constant Lagrange multiplier , in 

Eq. (26), may then be taken = t 


1. 4 Erdmann- Weierstrass Corner Conditions 

Eqs. (11) and (19) hold at any point on the extremal arc E(i,f) . 
Now, Eq. (11) indicates that at a corner, even if // is discontinuous, 
the canonical variables are continuous, i. e. , oj; - f V . ) . The 

(— ) and (+) signs indicate limiting values approaching the corner from the 
left and from the right side. Similar reasoning can be applied to Eq. (19). 
Thus, the following continuity conditions may be derived: 

W-A/ . <7'-^ 


(28) 
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Consequently, from the preceding continuity conditions it follows that at 
the corners 

/ / / ) ~ / , . \ + 

(29) 


("*) - ( h j>) 


if we assume that all of the state variables are continuous along the extremal 
(e. g. , it is assumed that at corners only the mass flow may be discontinuous 
and that there is no staging). Conditions (13) and (29) lead us to the con- 
clusion that at any corner point f i. e. , point at which flft) 

discontinuous ). - 0 . Based on this conclusion , the 

function H Q ft) * ft) ) ft) ] ~ //~ ('?) will be 

y 9 T i 

called the switching function since for 


is 


a) H~(T) i O , ft 


max. 


H s (T) = O , j& - var. or CoT^NEH 


C) H~(T)<o . ft-/S m , n . 

thus indicating when the mode of control shifts from one to another form 


along the extremal. 


ill 


2. GENERALIZED EXPRESSION FOR THE OPTIMUM BURNING PROGRAM 
ALONG THE /3 - VAR. SUB-ARC. 

As seen from Eqs. (13a), (18) and (20), along the - var. sab-arc, 


the following equations must hold 


H~ = - a 3 = O 

^ ■m ■ 



(30) 


(31) 


Thus, from total differentiation of Eq. (31) and accounting for Eqs. (1) to 


(3) and (14) to (16), it may be obtained 


z*. {*1L + 

m \ 'd Z 



Eqs. (31) and (32) lead to 


+ 3) - 5- + Cm =0 


( 32 ) 


(33) 


After total differentiation of Eq. (33) and using the preceding relations, 
we obtain 



(h + *)- s -*](4r % -A )+f‘[*( S 2 i+%)+ 3 ](V-) 
+.Af i ( S aS* S t)- h ]* - f C -A)/ = 0 


( 34 ) 
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Introducing the function 


(j> = x(S z +5)- S-m 


(35) 


it follows from Eq. (33) 



(36) 


Replacing Eq. (36) in Eq. (34), rearranging and eliminating C! , we see 
that the g eneralized optimum control program or g eneralized burnin g 
prog ram is explicitly given by 



The control program obtained in Eq. (37) is valid for any atmospheri c 
scheme adopted as well as any arbitrary aerodynamic character istics 
assumed, As will be shown, the control program y3 -= yS fej h j 272/ 
given in Eq. (37) is applicable for given-time or for free-time problems 
(i. e. , for e# o or C — 0 ). In fact, for problems where C = Oj 
the compatibility condition of Eqs. (31) and (32) is 



% 

'Tn 



(38) 
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which implies Cf> = o . Consequently, for free-time problems 
(C=o) , the optimum burning program along the y<3 - var. sub-arc 

may be readily obtained from the generalized expression in Eq. (37) as 


ji = MV SJ(S + *)- Szfcfyj; l %)- %) ,39, 

5 l (B+*V g + \ z ) 

Eq. (39) may be easily verified by taking the total differential 


jL y(*j,St)-.o 


(40) 
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3. TYPICAL VARIATIONAL PROBLEMS AND BOUNDARY CONDITIONS 
During the vertical ascent of a sounding rocket or during the phase of 
vertical ascent of a rocket booster prior to the initiation of the pitch-over 
maneuver, different optimal requrements may be imposed. For example, 
the rocket may be required to climb up to a certain final altitude with 
minimum propellant expenditure as it carries scientific equipment of 
maximum possible weight, or for a given propellant weight a maximum 
final altitude may be desired. In other cases, it may be desireable to 
attain maximum energy per unit mass when the given amount of propellant 
is expended (end of the vertical climb phase of a rocket booster) before 
starting the pitch-over maneuver. This may be of interest when large 
paylaods are put into orbit. At any rate, due to the zero- length launching 
conditions (vertical launching from rest) of large rockets, for ballistic or 
orbital missions, a necessary first phase of vertical climb, to which the 
problem treated in this paper is applicable, will always be required. During 
this phase it may be of interest to optimize a certain functional which is 
specified according to the optimal mission criteria chosen. Since this 
is a matter of the particular case considered and of the optimality criteria 
decided upon, we will only consider in the following work some examples 
of optimal problems that may be derived from the generalized formulation 

& * 6 ( Z J , z z , ,l F ,™ F , ? F )= min. 
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and to which the results previously obtained apply. The object here is to 
show how additional natural end- conditions may be derived depending on 
the minimal problem proposed, and boundary conditions imposed. 

3. 1 Maximum Altitude, Free-Time, Given Propellant Mass 

For this case the following boundary conditions will be assumed given 

V, - - «, - O s z x - C x = O 

- Zj - - O V 5 m - d, .O < 4 » 

Y 3 = / = o Yg * in F - e F - o 

^ j > ^ i > > dp y € p — Const. 

The function to be minimized is now 

G — hj - hp = bj - hp = frifi. ( 42 ) 

Consequently, from the Transver sality Condition and the boundary conditions 
(41) it may be found that 


( ^2 ~ dhp — O 9 dh p d’ 0 » • “ / 

CdZp = 0 y dtpjLO C-O-t 


( 43 ) 
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3.2 Minimum Fuel Consumption, Free- Time, Given Final Altitude 

Assume, in this case, that the following boundary conditions are given 

V s - «, = 0 fj-Cx-O 

= >>X ~ b I = O V 5 = - d F = O ( 44 ) 

y 3 = - 1 = o a ; ,b I , c z , d F = cW. 

The function to be minimized is 

G = 272 J - 772 F - 1 - 772 F = 7727/1. (45) 

Thus, from the Transver sality Condition, and conditions (44) and (45), the 
following natural boundary conditions follow 

di p = 0 , dz F + 0 ju = O 

~ t) dm F = O , dZ F ft O yU 3 - / <*<>> 

C dZ F = O , dZ F ft O 


C = 0= Const 
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Consequently, at the final point ~F it may be obtained that 



(47) 


Eq. (47) shows that for ^ O y JU 2 vanish at the final point if the 

arrival at this point is performed through a coasting sub-arc. For given 
time problems C 4^ 0 and then 



(48) 


3. 3 Maximum Final Energy Per Unit Mass, Free-Time, Given Fuel 
Consumption 

The following boundary conditions are now assumed 


^4 — z i O' i o 




'K • ~ *>! - o 


'K s ~d r -0 


(49) 


* x -i-o 


a z , c r ,d F = 


The function to be minimized is 


G = %)] = const. - d F (: + h 


= frtjn. 



118 


which assuming a normal extremal is equivalent to the optimal problem 

)win. ( 5 °) 




Thus, the following natural boundary conditions may be obtained 

( / i i fr ~ ^ Z F = O 9 ^ = ^F 


( ~ 0 d^F = 0 > dh F =£ 0 . * 


* / 


C dZ F - O , dZ F ±0 C- O - Const 


(51) 


Consequently, it follows that, at the terminal point T* 

(l-± 

1 J3 


\ - 


2 

/ rn 


If the final time were given, then C — Const. ^ O 


and 


( 52 ) 


A 


= rL fc - *£ 


z 3 




- ) -h 

J ™ , 


( 53 ) 


Eqs. (52) and (5 3) indicate that no arrival at the end-point V can be 
performed by a coasting sub-arc (Z = °) otherwise — 


=r OO 
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4 . THE AERODYNAMIC DRAG FUNCTION AND CONDITIONS ALONG THE 
y 3 -VAR, SUB-ARC 

In the analysis developed in previous paragraphs an aerodynamic drag 
function of the general functional form = ^ ) » has been 

assumed. On this basis a generalized expression for the optimum burning 
program has been derived. The aerodynamic drag function wUl now be 
given two specific forms in order to analyze particular applications of previous 
results. These forms will be used later for the numerical solution of 
some examples. In this analysis we will refer to free-time problems 

(c=o) ■ 

For the sake of discussion arbitrary aerodynamic drag functions of 
the forms 

a ) 3- Const 

(54) 

i) 3 -£ t z* = 3(z,h) , 

will be considered. These expressions, however, have some engineering 
meaning. For example, the form (54a) corresponds to the case of flight 
in an atmosphere of constant density and has been applied in previous 
investigations (i. e. , Refs. 6 , 7, 1 1). The form in (54b) corresponds to 
assuming the hypothetical case of a rocket having constant zero-lift drag 
coefficient and flying in an exponential atmosphere. 
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4. 1 Case 


£> = 6 f z 2 - T> fz) 


In this case, Eq. (38) leads to 


in = $ ( /+ z) = m (z) 


( 55 ) 


From Eq. (39), the optimum burning program or control program is 

% ^ + Z + 3 Z ) ~ \ 

ft = D -p f - = ft) 

2 + Z (4 + z) y 


(56) 


Finally, from Eqs. (l) f (55) and (56), it may be derived that the acceleration 


along the -var. sub-arc is given by 

>_ Z 3 + 3z 2 + 2 z 

ZZ 


Z + 5 Z + 6 z +2 


= Z’ft) 


(57) 


Eq. (57) shows that along the y3 -var. sub-arc the vehicle decelerates, 
which is consistent with what is implied in Eq. (55) since On must 
decrease. . 

In this case the numerical solution of practical applications is simple 
due to the two-dimensional character of the closed-form expressions (55) 
to (57). By the same token, and as will be shown later in the examples, 
the determination of corner points is straight forward. 

— 2 - ^2 h ~ ~ 

4. 2 Case T> = jj, g e _ = 

From Eq. (38) is now obtained 


in = D(1+z) = m(zft) 


(58) 
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The optimum burning program or control program is obtained from Eq. (39) 


as 


~ 5 (z+ s *)(?+*■)+ 4 z 2 ( > + zf 


2 + Z (4 + z) 


=fi(z,h) 


(59) 


The optimum acceleration along the a -var. sub-arc is obtained from 
Eqs. (1), (58) and (59) as 


( 2 + Z (4 + Z)J (l + z) 


~k z Z 3 + ( & 2 ~ 0 Z 2 - 2 z 

Z 2 + 4 z +2 




(60) 


In this case, the solution of numerical applications is more complex than 
in the previous case due to the tri-dimensional character of Eqs. (58) and 
(59). It may be readily seen that for , Eqs. (58), (59) and (60) 

reduce to Eqs. (55), (56) and (57) respectively. As will be shown in the 
examples, the determination of corner points can now be done using the 
corner line. This technique will be discussed in more detail later. Along 
the J3 -var. sub-arc the equations of motion and Eq. (59) require the 
use of numerical methods of integration. 
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5. EXAMPLES - MAXIMUM ALTITUDE FOR GIVEN FUEL CONSUMPTION 
AND FREE FINAL TIME 

The object of this paragraph is to show the application of the previous 
theory to the numerical solution of a given boundary- value problem. In 
particular, for the case where D = S(x,h) Special emphasis 

is placed on the determination of corner points and on the integration of the 
admissible set of variable Lagrange multipliers, or adjoint variables. 

For the problem proposed, the following boundary conditions will be 
assumed: 


V, 2 = o 




% - X - o 


X s z F -o 


( 61 ) 


V 5 = I -1=0 


% = m p - 0.4 = O 


The functional to be extremized is written 


G = hj -h F = h F — tnin. 


(62) 


As indicated in sub-paragraph 3.1, for this problem the following natural 
boundary conditions are obtained 



C — 0 — const. 


For the first example, the following values will be taken; 


D = k f z z exp (- & z h ) 


A 


= 2 


tn ax. 


, *,-5 , H.87 

’ / ml „. = ° 


(63) 


(64) 
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As indicated in paragraph 1, £ Eq. (13)J, three types of sub-arcs are 
admissible; /* max . > " VOn and ' The sequence 

of these sub-arcs will be discussed in the following. Due to the launching 
conditions in (61), obviously no = 0 sub-arc may be started at the 

initial point I. Also, if the initial values in (61) are replaced in the function 
(35), we get 

f t - ,b5) 


from which we concluded that no ^3 — VGP. sub-arc may be started at I. 

In fact, this also can be readily verified from Eq, (60) which gives the 
acceleration Z along the — VQP. sub-arc for different values of the 
parameter a z This is shown in Fig. 2. 

Thus, at point I the only sub-arc that may be started is ^ fr)ax ~ ^ ~ Cons ^ 
(see Fig. 3). Starting at I with a /3 sub-arc, our next problem 

is how to determine the position of an eventual corner point. For that, we 
will make use of the corner line . Integrating numerically the set of equations 
(1) to (3) with y3 = / ft mn , the state variables %(Z) , h (Z) , 'ftl(Z) 

are obtained. The integration was performed using an I. B. M. 1620 computer 
and applying the Runge-Kutta step-integration method. Replacing the functions 

/v /V 1 

z(z ) and h(Z ) obtained along the J3 max sub_arc in the equation 


m = 2 fS 2 + D) - T) = to [z(z) ,h(7)J 


( 66 ) 


the resulting function *77? may now be plotted in the 2 ^ -plane. The 
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line so obtained in the ( fr ? y Z ) -pla ne, will be called corner line (this line 
must not be confused with loci of corners), and is shown in Fig. 4 identified 
by the letter £ . Since at corners the state variables Z(Z') i h ( ) and 
K(V) are assumed all continuous (see Fig. 5), a corner point may only 


exist at the intersection of the corner line X and the sub-arc. 

y max 

Several examples are shown in Fig. 4, in which boosters delivering different 
thrust levels, i. e. , /3_ =°o, 5, 2 and 1.5 have been assumed. The 

y max. 

line A ^ A 2 A A 5 in Fig. 4 is the loci of corners mentioned before. In the 

particular case ^ ere considered, the corner point is identified 

in Fig. 4 with the letter A . The corner line ~t was extended past the 
point A in order to verify the non-existence of other possible corners. 

The line 2^ , within the range of interest, is monatomically increasing and 
therefore only one corner point appears admissible. In the numerical 
calculation of the corner point A , an automatic stopping condition 
included in the computer program as 


was 


| < ‘A.* r - ** ex F (~ ) ■(*+') I- c 


(67) 


where € > O is the precision with which the corner point A (Fi g. 3) is 
desired and , /)(Z) are the values along the J^rnax ~ ^ sub-arc. 

At the point /A , the transition to the y3-var. sub -arc was made, as shown 
in Fig. 3. Along the sub-arc, the set of Eqs. (1) to (3) with ^3 

replaced from Eqs. (39) or (59) was integrated up to the point £> (Fi g- 3) 
where On = = 0.4 . No departure from the ^3-VGr. sub-arc may 
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exist between A and -B , since this sub-arc is interior to the region of 
admissible displacement determined by the combination 
This combination satisfies the prescribed boundary conditions and 
within this region, both-sided control variations <f/3 = 0 leading to 

= 0 , ( t e., j3- vor.), are admissible. From point 3 the vehicle 
coasts to F (Fig. 3) to meet the boundary condition = 0 ■ It is easy to 
see from Fig. 3 that depending on the final values of 2 ' and/or tft p 

either one of the last two sub-arcs or both of them may not exist. 

A graphical representation of the extremal arc in the space of state 
variables (*,!> > ™ j is shown in Fig. 5. In order to verify the extremal 
properties of the arc IABF , position of the corner points and sequence of 
sub-arcs, the adjoint variables or Lagrange multipliers may be integrated 
backwards. Since at the final point 


[a ($•')] - A ,', 


( 68 ) 


and Z — 0 > then it follows that M = 0 ■ 

F / ’/r 


Now, Eqs. (14) and (15) 

may be integrated backwards (F-3) , yielding “ d A‘ (t) 

along the J3 — 0 sub-arc (F3 sub-arc in Fig. 3). At the point ^£> , the 
multiplier J^3 * s ^ oun< ^ using Eqs. (13a) and(l8). Now, Eq. (16) may be 


integrated forward again to calculate 

r 


A/ rj 7 Afif/~ s ) dz + A 


<69) 


B 
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along the sub-arc. Backwards integration of the system of 

Eqs. (14), (15) and (16), from the point 2$ > permits the calculation of 

all the adjoint variables associated with the extremal lAUF . The result 
of the numerical integration performed for the present example is shown in 
Fig. 6, where the multipliers and are 


shown. 


Note that at the corners A and _3 , the function A^ is 

discontinuous due to the discontinuity 4 ® in the control variable. Finally, 
the function may be calculated along the extremal and is shown in 

Fig. 7. Note that, as indicated before, for hj~ = 0 , y3 '> f° r 

= 0 , j$ « var. and for W~SO, / = /*,.„. A, corners A and 

3 , H s = 0 ■ 

Finally, the same boundary-value variational problem was solved using 
a drag function X) — Z , — S = Const In this case there was 


no problem in the determination of the corner point A since the Ji-Vor. sub- 
in the (m y z) -plane may be drawn a priori. In fact, from Eq. (55) 


arc 


On = 5 z 2 ( 1 + zj 


(70) 


along the 


/ 


-VQn sub-arc. The solution of the problem in this case is 


shown in Fig. 8. 
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LIST OF CAPTIONS 

Fig. 1 Hamiltonian H in terms of the bounded control variable yS 

Maximality Condition. 

Fig. 2 Acceleration along the J3~var. sub-arc. 

pig 3 Broken extremal solution of the proposed boundary-value 

problem. Case D= £, Z 2 <?•*/> (- & z h ) - 
Pip. 4 The /3 sub-arc for different thrust-levels, the 

o' J lYIQJC. 

corner line and determination of the loci of corners. 

Fig. 5 Broken extremal solution in the space of state variables 

(z, h, m) • 

Fig. 6 Lagrange multipliers associated with the extremal arc. 

Fig. 7 The switching function along the extremal solution, 

pig. 8 Broken extremal solution of the proposed boundary- value 

problem. Case D — ^ 2 
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Figure 1 
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Figure 6 
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Figure 7 
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NOTATION 


Elements of elliptical orbits: eccen- 

tricity, inclination, semi-latus rectum, 
period, perigee angle, respectively. 
Subscripts 1, 2, 3 refer to initial, _ 
target, and transfer orbits, respectively. 


Time intervals on orbits corresponding to 

ft 2 ’ and ’ res P ectivel y 


Relative time of nodal passage (positive 
when target passes before ferry) 


Angle between departure and arrival on 
transfer orbit 


Angle between departure point and node or 
reference 


Angle between arrival point and node or 
reference 


True anomaly of arrival point in transfer 
orbit 
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RENDEZVOUS POSSIBILITIES WITH THE IMPULSE 


OF OPTIMUM T./O-IHPULSE TRANSFER 


By 

D . F . Bender 



Summary 


Two-impulse optimum orbital transfer will lead to rendezvous for 
only a single value of the relative positions of the ferry and target 
at the beginning of the maneuver. It will be shown that, by the sim- 
ple expedient of splitting either the first or second impulse of this 
optimum transfer between any two orbits and holding for one revolution 
in the intermediate transfer orbit so obtained, rendezvous is possible 
over an extended range of relative phases . A technique for generating 
the required data including excess thrust tolerances for any orbit pair 
will be explained and sample data presented. 


I . INTRODUCTION 


In a large number of cases optimized two-impulse transfer between 
two orbits around an attracting center requires less total impulse 
than any other type. If a rendezvous at the conclusion of such an or- 
bit change is required, then it is clear that in general only one value 
of the relative positions of the two objects at the beginning will be 
allowable. It is desirable to discover the widths of such minima so 



14 



PROJECTION ON UNIT SPHERE 



FIG. L TRANSFER GEOMETRY 
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that error sensitivities can be indicated. In addition, it will be shown 
that by permitting three impulses the minimum of the total impulse versus 
relative phase for any pair of orbits can be made to have a horizontal 
base whose width is at least the difference in periods of the initial and 
final orbits. The technique is to split one of the two impulses into two 
parts used exactly one revolution apart. 

The geometry of the transfer is indicated in Figure 1, in which the 
orbit planes are projected onto a unit sphere. If the orbit planes of 
the passive target satellite (terminal or 2) and the ferry vehicle (ini- 
tial or l) are inclined as indicated, the line of nodes is taken as the 
reference direction with the ferry ascending. If not, one uses ary di- 
rection in the plane, usually the perigee of one of the orbits. The 
relative phase of the two is given by indicating the position of the tar- 
get (2) when the ferry (l) crosses the reference direction line. In this 
discussion the time interval t is used, that is, the target is at a posi- 
tion corresponding to the time r past the reference line when the ferry 
crosses it. 

Two-impulse transfer has been expensively studied as a perusal of 
the aerospace and astronautics journals will indicate. For this discus- 
sion it is necessary to have a computer program which is able to survey 
and optimize on total impulse so as to select from all the possible 
transfers those requiring the least fuel. Such a program has been devel- 
oped by Kerfoot and DesJardins^ and further improved by Me Cue 2 of the 
Space Sciences Laboratory of S&ID, NAA. 

The formulation of the two-impulse transfer problem by DesJardins 
and Kerfoot is to express the total impulse as a function of three 
angles ($n, j $->'), or (^). The angle is the variable which selects 
a particular transfer orbit between the departure point (^ ) and the ar- 
rival point (ftp) . Hie variable used is the true anomaly of the arrival 
point in the transfer orbit considering AO < v . In all cases both a 
short (AS < tt ) and a long (AS > v) transfer are considered and the 
better one selected. (The discarded transfer corresponds to motion in 
the opposite sense on the transfer orbit . ) 

The f ull range of phasing possibilities is encompassed by 0 < r < P 2 
(xdiere ?2 is the period in the second orbit), since any value of r out- 
side this range may be treated modulo P 2 * Consider next the ensuing revo- 
lution of the ferry. The value of r has decreased by the difference 
P 2 - Pq, and it is clear that on successive orbits of the ferry the value 
of r will continue to step by this difference. 


If tq, t 2 , and t-^ are the traverse times associated with the true 
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anomaly intervals ^5q, and A0 then a necessary and sufficient condition 

for rendezvous is tq + to = t£ - T or T = t^ - tq - tj. Since tq, t 2 , 
and to are functions of it follows that T is also a function of 
Thus, 'rinding optimum two-impulse rendezvous trajectories is a matter of 
minimizing l(j0) under the constraint = constant. It is easy to see 
that this consists of finding points in jfi space at which the surfaces 
I = constant and T = constant are parallel. Analytical expressions for 
all the required derivatives can be obtained easily3 but are not given 
here. They were programmed and used in the searching technique for the 
time constrained optima. 


II. THE IMPULSE SPLITTING TECHNIQUE 

The relative phase, t, between two vehicles cnanges by the diiference 
in periods each time the initial vehicle crosses the reference axis. If 
a vehicle holds in the initial orbit, the decrease in r can only be 
P 2 - Pq. If, on the other hand, a hold is possible in some intermediate 
orbit of period P 1 , the effective value of T will decrease by Po — P • 

By splitting either Iq or I 2 into two parts, the second of which is used 
one revolution after the first, any period P between Pq and P 2 may be 
attained. Thus any value of r lying between an optimum r 0 and r D + (Pq - p l' 
is effectively reduced to r 0 and made accessible for rendezvous. This 
technique is similar to the looping methods described by Silber and 
Hoelker.^ Two immediate consequences are evident: 

(l) An upper bound is provided on the number of i-raiting periods 
in the initial orbit before an optimum rendezvous maneuver 
may be performed. The bound, n^^, is given by 


n max - 



(2) If as many as + 1 revolutions are permissible, it is 

then true that athree -impulse rendezvous maneuver can be 
made which requires no more fuel than the optimum two- 
impulse orbital transfer. 

For the case when Pq — P^ — p 2s P* H^st hie in the range Pq to Pq» 



and thus the curves of constrained optima are translated from the positions 
indicated in Figures 3 and 6 to points (P2 - Pq) greater and the minimum 
value of impulse may be achieved anywhere in the interval. If the trans- 
fer orbit period happens to lie outside the range Pq to P 2 , greater ranges 
of P 1 and hence of r are accessible. If P 3 < Pq the range for r extends 
to P 2 - P 3 above that of the constrained optimum. If P^ > P£, the range 
is from P 3 - T '2 below to V 2 - Pq above. Both of these cases occur in 
Figure 6 . 

Properties of the intermediate transfer orbit depend upon the geometry 
of the transfer problem and the fraction of impulse used in the phasing 
maneuver. Velocity components may be determined; then angular momentum 
and energy may be obtained and used to calculate the elements of the inter- 
mediate orbit. 

Up to this point discussion has been limited to multiple holds in the 
initial orbit. In ary case where P3 lies outside the range Pq to ? 2 , 
fewer total holds may become possible if more than one is taken in an 
intermediate transfer orbit. 


III. TYPICAL IIUKKRICAL RESULTS 

numerical results using this technique are shown for two different 
pairs of orbits in two sets of figures: Figures 2, 3, and 4 for a pair 
of circles inclined at 3-5 and Figures 5, 6 , and S for a pair of inclined 
and asymmetrically oriented ellipses. 

The two circular orbits of Figures 2, 3, and 4 have radii of 4070 
miles and 4270 miles respectively. The optimum impulse versus departure 
point is shown in Figure 2 where the effects of the inclination are evi- 
dent. Since there is no orbital distinction between the two nodes, only 
one optimum needs to be explored in the search for impulse versus relative 
time ( r ), Figure 3. The stepping of r can be considered to be with steps 
of (P 2 - Pq)/ 2 every half revolution of the ferry. 

These two circular orbits are close together in period and illustrate 
the need for proper phasing if rendezvous is to be accomplished in a few 
revolutions, since even the extention of the minimum of Figure 3 by 
P 2 - Pp to the left covers only a small fraction of the total period, P 2 . 
In Figure 4 the impulse splitting possibilities are illustrated for the 
optimum and for three points on the curve of Figure 3 requiring slightly 
greater total impulse. The excess is indicated as a percentage on each 
line. If the line slopes upward to the left the first impulse is to be 




DEPARTURE POINT, 9„ DEGREES 

FIG. 2. IMPULSE vs. DEPARTURE POINT 
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split, but if it slopes upward to the right the second is to be split. 

The whole range of r which has to be covered is Pp anc * as shown in Figure 5 
this may require as many as fourteen orbit holds before initiating the 
rendezvous portion if the worst possible phase relation should develop as 
an initial condition. For the optimum case, the total impulse and the 
values of impulse at departure and arrival are indicated. 

The second pair of orbits are two inclined ellipses with elements: 


p^ = 5,000 mi. 

ei — .2 

1 ~ 

-90° 

il = 5 

Pp = 6,000 mi. 

e 2 = * 2 

2 = 

0 

O 

V 

h ~ 0 


The possibilities for rendezvous are much better. In the first place 
there are four optimum two-impulse transfers as shown in Figure 5* The 
two orbits would intersect if in a plane and actually pass quite close 
to one another at the ascending node. Each of the optimum is charac- 
terized by whether it is short or long (A 9 < it or > tt) and internal 
or external to the two ellipses. The latter distinction applies also 
to the periods, that is, internal has < Pq and external, Pjj > Pp. 

In Figure 6 the results of the time constrained search are indicated 
for each of the four ojPbina as well as the range of relative phase ( r ) 
accessible to rendezvous by splitting one of the impulses. For any point 
of Figure 6 a range of P 1 lying outside the interval Pp to Pq Is possible 
in two ways, that is, either the first or the second impulse may be split. 
The nature of the curves is illustrated in Figure 7. 

Finally, in Figure o the splitting fraction is plotted for each of 

the optima and for the six encircled points of Figure 6 except that 

only one of the impulse splits is indicated over the overlapping range. 

The percentage impulse required over the lowest is about 15$ • Again the 
total impulse for each optimum and the impulses at departure and arrival 
are Indicated, and which impulse is split can be seen by reference to 
Figure 7* For this case only a harrow range of r (l, 700 seconds to 
2,300 seconds) is not accessible to three-impulse rendezvous within the 
1 5$ limitation. For phases in this small region, a hold of one revolu- 
tion in the initial orbit is required. 
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NOTATIONS 


A 

A 

— E 
E 

E ' 

AE 

J 


J* 


J 

J 

J 


1 

2 


J** 
J ' 
l 

L 


origin of rotating coordinate system 
position vector from barycenter to center of the 
rotating system 

the position vector of A relative to the earth 

position vector of the earth relative to the bary- 
center at t = 0 

position vector of the earth relative to the bary- 
center, but rotated through an angle ccT 
E ' - E 

Hamiltonian (Jacobi integral) for the restricted 
problem 

difference between the restricted Hamiltonian 
and the two-fixed-center Hamiltonian 
the part of J independent of a, ft, and y 

the part of J that is a function of o', ft, and 7 

Hamiltonian equivalent to J* but written in terms 
of two-fixed -center coordinates and momenta 
time dependent part of J* 

Hamiltonian of two -fixed-center problem 
length of position vector from earth to moon 
position vector from earth to moon 


L 

L 


position vector of the moon relative to the earth 
in the rotating system 

velocity of moon with respect to the earth (O x L ) 



fl x L in the rotating system 

momentum canonically conjugate to 

momentum canonically conjugate to 

position vector relative to a point fixed in inertial 
space e.g. barycenter 
position vector relative to the earth 

position vector relative to the moon 



position vector relative to A in the rotating system 

position vector relative to the earth in the rotating 
system 

position vector relative to the moon in the rotating 
system 
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NOTATIONS (Cont'd) 


— E 
R 


— m 



a 

8 

y 

6 


u 

n 


CO 

grad y 


= position vector from barycenter to earth 

= position vector from barycenter to moon 

= position vector relative to A in the rotating system 
for the two -fixed-center problem 
= position vector relative to point at A 

= length of position vector relative to earth 

= length of position vector relative to moon 

= a specific period of time 
= time variable 

= constant coefficient of L. in composition of A 

= constant coefficient of O in composition of A 

= constant coefficient of Ij in composition of A 

= the angle of rotation of the coordinate system 
about the barycenter after a time T 
= gravitational constant of the earth 
= gravitational constant of the moon 
= angular velocity vector of the moon about the earth 
= the magnitude of angular velocity Q 
= gradient with respect to the components of V taken 
as coordinates 


Subscripts 

B 

0 


vector relative to the barycenter 
initial value 


Superscript 

dot over quantity = first total time derivative 

2 dots over quantity = second total time derivative 
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Approximation of the Restricted Problem 
by the Two-Fixed Center Problem 

By Mary Payne 


SUMMARY 

. te* 

In this report, a perturbation theory of the two-fixed-center problem 
leading to an approximation for the restricted-three-body problem is developed. 
It makes use of a generalization of the method developed at MSFC by Schulz- 
Arenstorff, Davidson, and Sperling. ( 1 ) The derivations are carried out in a 
coordinate system rotating about an accelerated origin, and the generalization 
consists of the selection of this origin in such a way as to minimize the effects 
of the non-integrable terms in the perturbation equations. The results of some 
numerical calculations are presented. 

INTRODUCTION 

The equations of motion for a vehicle moving in the gravitational fields 
of the earth and moon are: 



( 1 ) 
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where R^, R£, and R are the position vectors of the vehicle referred to the 

earth, the moon, and a point fixed in inertial space, respectively. Lower case 
letters denote the magnitude of the corresponding vectors. In this report it 
will be assumed that the earth and moon are moving in circles, under their 
mutual gravitational attraction, about their common center of mass. This prob- 
lem is the restricted three-body problem, and the fixed point may be taken to 
be the center of mass of the earth and the moon. An approximation to the solu- 
tion of the restricted problem will be sought in terms of the known solution!* 1 ) to 
the Euler problem of two fixed centers of gravitation. The method will, in many 
respects, follow closely that developed by Schulz-Arenstorff, Davidson, and 
Sperling. In their procedure, the problem is transformed to a coordinate 
system rotating about the center of mass. In this rotating system, the Euler 
problem is taken as the basis of a perturbation theory. Using the initial con- 
ditions of the Euler problem as a set of canonical variables, it is shown thati > 


R = + grad J* 

u -0 

and . (2) 

-0 * - grad R 0 J*. 

where R Q is the initial position vector in the rotating system, P Q is the momen- 
tum vector conjugate to R Q , and J* is the difference between the Hamiltonian for 

the restricted problem (Jacobi integral) and that for the Euler problem, and is 
given by 

J* = fl • Rj x Pj + J** . ( 3 ) 

The solution of the restricted problem is given in terms of an osculating two- 
fixed center problem with varying initial conditions. If J** were zero, the 

equations for R Q and P Q could be integrated directly. In the Schulz-Arenstorff 

theory, J** does not vanish and, in fact, contributes appreciably to the vari- 
ation of Rq and Pq if the time interval over which the integration extends is too 

large, or if either the earth or the moon are approached closely by the vehicle 
during this time interval. 

It is the purpose of this report to show that the effect of J** can be reduced 
by selecting an origin for the rotating system other than the center of mass of 
the earth and moon. In the course of this development the details of the Schulz- 
Arenstorff method will be given, and the coordinates for a center of rotation will 
be determined so that J** and its first time derivative vanish initially. 



PRELIMINARY CONSIDERATIONS 


Since the two-fixed center problem will be used as the basis of a pertur- 
bation theory, it is necessary that the earth and the moon be fixed in the ro- 
tating coordinate system. This implies that the origin of this rotating system 
must be fixed relative to the earth and the moon. The most general of such 
points will rotate about the barycenter with the angular velocity of the earth and 
the moon. The radius vector from the barycenter to the origin of the rotating 
system can be expressed as 

A = O' L + 0Q, + yh , (4) 

where L and L are the position and velocity vectors, respectively, of the moon 
relative to the earth in a non-rotating coordinate system, and Q, is the angular 

velocity of the moon about the earth. From the definition of L and L it is ap- 
parent that both vectors are known functions of time. Furthermore, L and L 
are constant vectors in the rotating system and ^2 is constant in both the 
inertial frame and the rotating system. Thus, the requirement that the point 
A be fixed relative to the earth and the moon implies that a, 0, and y are 
numerical constants. The constant 0 may be chosen arbitrarily, for the point 
A is used to determine an axis of rotation oriented in the C2 direction, and all 
points with the same a and 7 will lie on the same axis independently of 0. Thus, 
0 may be taken as zero without loss of generality, and it will no longer appear 
in the formulation. Referring to Figure 1, it is seen that R, R^, _R 2 , L, and 

R^, the position vector of the vehicle relative to A, satisfy the following re- 
lations : 


— E 




*1 


—A 


—A 


- -S-' 

M + M 


*2 =± 


L 

L 


= R, + - A = R - (a+ — f 4 -*) L - y L 

—i — e — —i v n+ /Liy - — 

= R, + R. . - A = R 0 - (a — L - y L 
—2 — M — —2 V \l+ [l s — — 

+ R . - R ^ + o; Lj + 7 L 


(5) 

( 6 ) 

(7) 

( 8 ) 

(9) 


R 


A 


(10) 
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First, it is necessary to eliminate R from Eq. (1) and obtain the equations of 
motion in terms of R^, R^ and R 2 - To do this, one may differentiate Eq. (10) 

twice with respect to time: 

• • «• •• * * • 

R = R A + a L + y _L . (11) 

Now, the condition that the earth and moon move in circles under their mutual 
gravitational attraction means that 


L = fl x L 


and 


L = x L 


(M + m' ) “3 

l 


(12) 


Differentiation of Eq. (12) (with i = 0, as L has constant magnitude), enables 
us to write Eq. (11) as 


n - r - H tJ L 
— —A .3 


a L + 7 


-] ’ 


(13) 


and the equations of motion (1) become 


*A = ' ^ “3 ' ** 


Si 


R 0 , / 

/ z2l , M + M. 


(«l + y l ) . 


(14) 


It should be noted that, at this stage, the coordinate system associated with A 
is an accelerated system since the origin has uniform circular motion. It is, 
however, not a rotating system yet - that is, the coordinate axes remain parallel 
to the inertial axes at the barycenter. 

The next step is to transform to rotating coordinates about A. The vectors 
in this system will be denoted by bars, and the equations of motion become 

R^ = - (it — - p,' — ^ 2 *^ ( 0 ! L + "y L) ~ 2 O x R^ - Q x x R^^) . (15) 
^ T 1 V 2 1 

It should be noted that, in this rotating coordinate system, the earth and the moon 
are fixed, with position vector Ij of the moon relative to the earth as a constant 
vector. The vector L does not represent the velocity of the moon (which is zero), 
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but is a vector mutually perpendicular to L and Q , and satisfying Eq. (16) with 
bars over the vectors. As the rotating system has angular velocity O , it follows 

of course, that Q and £2 are identical. 

A constant of motion for the problem in the rotating system may now be 

obtained by dotting Eq. (15) with and noting that the earth and the moon are 
fixed in this system, so that 




-i 1*2 ' 


(16) 


Thus, 


5 a ' 5 a 


— 2 

, .r; 

d f —A 
dt V 


¥■)• - 


Mi 


-i‘ 5 X 


, - 2 ’ -2 


(17) 


lU j u ~ (a R a - l + yR A - L)+(pxR A )- (qxrJ 


+ ,t V L '(“5Ai + ’'5A-i) + *(0 x B A y 


r l ‘2 l' 


as L and L are constant vectors. Denoting the constant of motion by J: 

2 2 
j = *5a - - lt ^ i - (“BA-S + ’'B A -i : )-4(nx| A ). 


(18) 


It may now be shown that, if the vector 


P A = R A + 0 x R a 


(19) 


is regarded as the momentum conjugate to R^, the integral J of the motion be- 
comes the Hamiltonian. To prove this, substitute for R . , using Eq. (19), in 
Eq. (18): 

2 2 

J=K£a-S x Ia) (20) 


= *?A 2 -f 2 -a-BA^A-^^SA-i^BA-i 
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If R and P . are conjugate vectors, Hamilton's equations, 
“A A 


R* = S r& d t) J = P A " 0 x R. , 


-A 


—A 


( 21 ) 


and 


-A 


- - s rad g J - *-| A + £> - a * £a + C° i + r i ) <22) 


must be satisfied. It is evident that Eq. (21) is identical with Eq. (19^ defining 
the relation between velocity R A and the momentum P A conjugate to R A - Now, 

it will be shown that Eq. (22) reduces to the equations of motion (15) in the ro- 
tating system. First, 


(23) 


But, 


S rad R £ = 2 Sradj; ^ . 

-A 1 r 1 -A 


> 


(24) 


hence, 


2 r l S rad R A r l ' * rad g A ^ 


■ « rad R A (^A - 5 e + a) 

2 

= gradg [ R A +2R a - (a - R e ) + (a - R e ) 

—A 

= 2R A + 2 (A - R e ) - 2R X , (25) 


so that, finally, 


and 


. -i 

8 rad R A r l = Ff 


(i 

gra % % - - T7 


(26) 
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Similarly, 


grad- = “ 

r 2 


4 Be 


3 5 


(27) 


so that Eq. (22) may be written as 


• 4 R-i u'R 0 , / , — v 

-A = 3" F " — x —A + U ,3 M C® — + y — ) 


(28) 


1 2 

Now, from Eq. (19), 


-A " -A + ^ X -A ’ 


(29) 


and use of this relation for _P^ and Eq. (19) for in Eq. (28) yields 

4^ u'R 2 


R a + Q x R a = 


3 3 11 x R a — (1 x x R a ^ 

r l r 2 


(30) 


“-y*-' (c,L + yi). 


Finally, if the Ox R A on the left is transposed to the right hand side of Eq. 

(30), it becomes identical with the equations of motion (15) in the rotating sys- 
tem. 


RELATION BETWEEN THE TWO-FIXED CENTER PROBLEM 
AND THE RESTRICTED PROBLEM 


A Hamiltonian, J, has now been obtained for the restricted problem in a 
rotating coordinate system with the origin at A: 



CL 


R . x 
—A 


-A - (“ -A i + ™A- i)- 


(31) 


with 


A = a. L + 7 L 


(32) 
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referred to the barycenter of earth and moon, and 


p a = B a + Q * B a • < 33) 

The development so far differs slightly from that of Schulz -Are ns tor ff, Davidson, 

and Sperling^ in two respects: it has been carried out in three dimensions in- 
stead of two, and the center of the rotating coordinate system is at A instead of 
the barycenter. Following their development, a solution of Eq. (31) in terms of 
the solution of the two-fixed center problem is now sought. For the two-fixed 
center problem, the Hamiltonian is given by. 


J' = 


/ 

_ -tL 
r * 

2 

(34) 

and the Hamilton equations 

are 


k 

= grad , J' 
-A 



and 

k 

= - grad- J ' 

-A 

Si' , £ 2 

- 3 ~ H 3 

r i r 9. 

(35) 


Denoting the solution of the two-fixed center problem by primes and that for the 
restricted problem without primes, the solution sought is to have the form 


a (a„. p 0 > 0 ' a' (— o «• foW’ 0 


£ ( b 0 . £ 0 ’ *) = S' (£<,<'>• 0 • 

Thus, the problem is reduced to finding the time dependence of the initial con- 
ditions in the solution of the two-fixed center problem that provide the solution 
of the restricted problem in the same functional form as that of the two-fixed 
center solution. 

The theorem, mentioned in the introduction, on the equations determining 
the time variation of the initial conditions will now be given a precise statement. 
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Theorem: If R(R Q , P Q » t) and P(R Q , P Q , t) constitute 

the solution of a problem with Hamiltonian J 0, P) while 
R 7 (Rq, P q , t) and P 7 (Rq , P Q , t) constitute the solution 

of a problem with Hamiltonian J 7 (R 7 , P 7 ) with 

R(5 0 «?0’ 0) = 5'(R 0 ’^)' °> = 5 0 

and (37) 

P(R 0 > P 0 * = p / ©o' P 0 * °> = 5 0 


then Eqs. (36) are satisfied with R„(t) and P Q (t), de- 
termined by the equations 

R Q (t) = grad p ^ J* (R Q , P Q , t ) 

and (38) 

P Q (t) = - grad^ J* (5 0 > P 0 . t) » 

where 


J <R 7 , P 7 ) = J (R 7 , P 7 ) - J 7 (R 7 , P^= J* (Rq, Pq, t) (39) 

Wherever R Q and P Q occur on the right hand side as a re- 
sult of the gradient operations, they are to be replaced by 
Rg(t) and P^(t), respectively. 

(2\ 

This theorem has been provenbyArenstorf ' in an unpublished note and will 
now be applied. 


To obtain the differential equations for R^(t) and P^(t), J must be 

written in terms of R 7 and P 7 . , associated with the two-fixed outer problem. 
That is, A A 


J =J (5a* ?a> " j '®a> ?a> 

= • SI x Ea K-Prli'Ei, 


(40) 


where J (R^, P ^ ) is obtained from Eq. (31) by replacing R^ and P^ by the 
corresponding primed quantities, and J 7 (R^, P^) is given by Eq. (34). 

It is now necessary to obtain J* by expressing J in terms of the initial 
conditions of the two-fixed center problem. This is very difficult to do ex- 
actly, as the solution(O) of the two -fixed center problem is given in terms of 
elliptic functions with the initial conditions entering not only in coefficients of 
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these functions but also in their moduli. Therefore, the solution of the two- 
fixed center problem is a transcendental function of the initial conditions. An 
approximate solution is, however, obtainable by expanding J as a power 


series in time: 

_ _ « t 2 
J = J (0) + J (0) t + J (0) — + . . . 

- t 2 

= J* (0) + J *(0)t + J* (0) ~2 + • • . 


(41) 


Using Eq. (40), the first time derivative of J is 

j = -n • r;x p;-o -g;xp^ -t^ (0 R A -L + >-s;' L). (42) 

Xj 

Now, Eq. (42) contains time derivatives of and P A> which may be 
eliminated by means of the Hamilton equations (35) for the two-fixed 
center problem: 

LtR,' juRo _ “ 

j= -Q* p , a x ^a "2 'l' x( 1 T )_ 

r l r 2 1 


(43) 


The first term in this equation vanishes. Evaluation of J and J at t 0 yields 
J(0) = J*(0) = -Q ■ Em * ? m - ^ IIo ’ < “ 1 + ? ± >• (44 > 

Xj 

and 


m'R' 


j ( o, - bm - - 9 -f A0 ~ f > -^^0 • ^ 

(45) 


r 10 r 20 1 


Setting 


h = -8 • 5ao x -ao 


(46) 


and 


J 2 ' 


-■^RaO* («L+yL) +t[a-R^ 0 x ( ^ + - 3 ) 


11R' 10 H'K 


20 


10 


20 

1 


- IL ^ L J + 


(47) 


so that 


J* — + J2 • 
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Application of the Arenstorf theorem, now yields 


and 



= +grad , J* 
-AO 


- o x r ; 0 + gradp^ J 2 



grad 5Ao J * = ‘- X2 A°" gra %o J 


(48) 


(49) 


as the differential equations for the variation of the two-fixed center initial con- 
ditions, which must be included in the two-fixed center solution in order that it 
may become the solution of the restricted problem. 

If Jg were zero, Eqs. (48) and (49) would integrate immediately. They 
would simply say that and rotate clockwise with angular velocity 0 . 

That is, in the rotating system the solution of the restricted problem at time 
T would be given by the solution of the two-fixed center problem at time T, 
with initial conditions obtained from those of the restricted problem by a 
clockwise rotation through f} T about the point A. For T=0, the restricted 
and two-fixed center problems have the same initial conditions and, hence, 
have exactly the same solution. 

Actually, of course, does not vanish, and it is here that the selection 
of the point A enters. Every term of J ^ involves either or P^ , which 

depend on the selection of the point A, so that this point should be selected so 
as to minimize the contribution of J ^ to the variation of the initial conditions. 

This could be done in various ways. Inasmuch as the position of the point A 
depends on the two parameters a and y , it is evident that only two conditions 
can be imposed on the selection of A. Several such conditions suggest them- 
selves immediately: 

( 1 ) 


( 2 ) 


( 3 ) 


The first method has the disadvantage that the validity of the approximation 
would deteriorate with time, and there is no obvious way of estimating the 
duration of validity. The other two methods have the disadvantage that, if the 
time interval specified is too long, the approximation would not be valid, even 
initially, and again, a criterion for "too long" is missing. It was, therefore, 
decided to try the first method, which would give some insight into the duration 
of validity, and might very well produce results of practical value. 


Determine a and y so that in the constant term and the 

coefficient of t vanish for the initial values of R' and p' . 

—AO -AO 

Determine a and y so that J ^ vanish for t=0, with initial 
values of Rj^ and P' , and also vanish at t=T, with the 
rotated values of R ' and P^q determined by at time T. 

Determine a and y §o that the square of J 2 is minimized 
over the time interval 0 to T, using either the initial values 
of R^ 0 and P' or their time dependent values determined 

by over the interval. 
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DETERMINATION OF a AND y 


In accordance with the conclusion of the last section, a and y are to be de 
termined by the equations 

R -(aL+yL) =0 < 50 > 

AO 


and 



M -10 

3 

r io 


R, 


+ M 


20 


20 


;) 


U + h' p 

- ,3 4o 


(aL + yL) =0, 


(51) 


so that the first two terms in the power series expansion of in Eq. (47) vanish. 
The primes have been omitted in Eqs. (50) and (51) because tne initial values of 
R' and p' , regarded as variable parameters for the restricted problem, are 
the°initial'values of the restricted problem by the Arenstorf theorem. < > Now, 

R and P depend on the selection of the point A, so that, for the determination 

of a and y from Eqs. '(50) and (51), they should be replaced by the position and 
momentum of the vehicle relative to some point independent of A. A particularly 
compact form is obtained for the equations of a and y by replacing P^q by 

and R by R 10 or R 2Q , as follows. First, since from Eq. (19) 


Eao=5 ao + 9*5 A o' 


for any point A fixed relative to earth and moon, it follows that 


2io = -io + ^ x —io • 


(52) 


(53) 


Therefore, since in the rotating system the velocity of the vehicle relative to the 
earth is the same as that relative to A (both are fixed points in the rotating system), 


-A0 




+ £1 x ( R aq - R 10 ) 


-Q * ( (“ + * 


i ) 


- ( a + tStp- t+y 


y+M - 


3 - 


L , 


(54) 
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on making use of Eqs. (8) and (12). Thus, the third term of Eq. (51) will be 
proportional to 


P A0 • (a L + y L ) = P 1Q • (OL + yL ) - ? j- , 


where the terms in ay have cancled out. 

Again using Eq. (8), the first term of Eq. (51) will involve 


(55) 


- -AO * -10 " " — X 


(a+ u / nT ) ^ + y ± _ 


R 


M +M' “ J -10 


and the second term will be proportional to 


fl- *ao x 


K 20 = -Ox[(a- 7 r f ir )L tyL ] • R 20 


= ->*"5 — i] ■ 

Xj 

so that Eq. (51) may now be written as follows: 

io'i ^^Sxo-S] 

r 10 1 


(56) 


(57) 


(58) 




or, collecting terms in a and y : 
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r : 


a L -10 ' ± ( 3^3 


a' LL ' , U +\l ' ™ T 

L ) + — 3 —10 — 


r 10 r 20 


U + U 


R • l- 1 

' r -io - 
- y\u 


+ v 


10 


/ -20 ' ~ 


20 


+ rL-p -l! 
+ SL -10 - J 


(58) 


u U ' s . T ‘ / 1 _ _J \ 

r 5in i 1 3 3 ' 


VT+ W -10 - ' 3 

r 10 20 


= 0 


where use has been made of the fact that 


-10 ' - -20 - 


(59) 


Using Eq. (8) once more, one obtains for Eq. (50): 

r 10 • (ffL + yL ) - [( a + - ^ t - )L +yL ]•[<* L+yL ] 

~ u' 2 2 ,2 

= R 10 • (a L +yL ) -a( a+-^j7)£ -y ■ * 


aV + a(Ri 0 - L 


(60) 


• L -T^TT £ 2 ) -y 2 ^ +y® 10 ‘ L ) = 0 


If Eqs. (58) and (60) are solved for a and y , a point A is determined so that the 
following procedure should give an approximation to the restricted problem valid 
for_a time interval whose length depends on the size of J* and the rate of variation 
of R' and P' . The procedure is carried out in the rotating system as follows: 


Modify the initial conditions of the restricted problem by a 
clockwise rotation through to T about the point A, and solve 
the two-fixed center_problem with these modified initial 
conditions. Then, R^(T) and P^(T), given by the two - 

fixed center problem, should match R^(T) given by the 

restricted problem with unmodified initial conditions. 
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APPLICATION OF THE METHOD 


In order to carry out a numerical test of the method, use was made of the 
Republic interplanetary trajectory program. The input for this program requires 
that initial conditions be given in a coordinate system with its origin at the earth 
and axes with fixed directions in space. The z-axis points towards the pole star, 
the x-axis points to the first point of Aries, and the y-axis is selected so that the 
system is orthogonal and right-handed. The output includes coordinates and ve- 
locities of the vehicle in this same system. An option is available which fixes the 
moon at any desired point on its orbit and computes a two-fixed center problem 
for this fixed position of the moon and given initial conditions. A set of initial 
conditions is available which yields a lunar trajectory (referred to, henceforth, 
as the base case) with a moving moon, starting near the earth, closely circling 
the moon and returning to the earth. Thus, to test the application one could 
modify the coordinates and velocities at various points on this base case and com- 
pute a two-fixed center problem from the modified conditions to obtain a com- 
parison, which should indicate the time intervals over which the approximation 
is useful for various portions of the trajectory. 

The modification of the initial conditions derived in the preceding sections 
was carried out in a rotating system, and it is now necessary to transform this 
modification for use in the coordinate system of the interplanetary program. To 
see how this may be done, suppose for the moment that the point A is at the bary- 
center, i.e. , a and y are both zero, and that the fixed and rotating systems are 
coincident at t = 0. It is evident, in this case, that the two-fixed center orbit 
obtained from the initial conditions, modified by a clockwise rotation through an 
angle 0 about the barycenter, is exactly the same relative to the earth and moon 
as if the initial conditions had been unmodified and the earth and moon had been 
rotated counterclockwise through 0 about the barycenter. Now, the angle 0 is 
U)T, where T is the time at which the comparison is to be made. Hence, if the 
earth, the moon, and the two-fixed center orbit, corresponding to the modified 
initial conditions, is rigidly rotated counterclockwise through wt, the earth and 
moon will coincide with their positions at time T in the fixed system, and the 
point corresponding to time T on the two-fixed center orbit is the one to be com- 
pared with the restricted problem carried out in the fixed system. Moreover, 
this counterclockwise rotation just transforms the two-fixed center problem, 
with modified initial conditions and earth and moon in initial position, into that 
with unmodified initial conditions and earth and moon in their T positions. There- 
fore, for O' and y both zero, the comparison can be made, using the interplan- 
etary program by fixing the moon in its T position and referring the unmodified 
initial conditions to the coordinate system centered at the earth at time T. This 
is indicated in Fig. 2, where the unprimed initial conditions are referred to the 
earth at t = 0, and the primed initial conditions refer to the earth at t = T. The 
initial conditions are fixed. 


173 


A comment on the relation between the momentum vector Pg, conjugate to 

Rg, and the velocity vector Rg, where B is used to indicate that the barycenter is 
the origin of the rotating system, in now in order. Recalling the definition of Pa 
in Eq. (19), it follows that 

Eb = 1b + 0 x 1 b • (61) 

and hence Pg is simply the velocity vector in the fixed system with its components 
referred to the instantaneous rotating axes. Since it has been assumed that the 
fixed and rotating systems are coincident at t = 0, it follows that 

P , (62) 

o o 

where Rj$ is in the fixed system (recall that bars denote rotating system). At 
time T7if the Pg vector is rotated through WT counterclockwise, it will become 
the Rg vector. But this is just the transformation that has been used to translate 
the two-fixed center approximation from the rotating to the fixed system. 

Thus, if the barycenter is the origin of the rotating system (i.e. , a = y = 0), 
the prescription for the approximation is the following: 


(1) Let 


(L(T)-L(O)) 

be the displacement of the earth in time T. 

(2) Set 

5' 10 =5io- ^ = %o + ^' 


and 


-10 -10 ’ 


since a translation of the origin will not affect the velocity. 


(63) 


(64) 

(65) 


(3) Fix the moon at L (T), that is in its position at time T relative to 
the earth. 

(4) Solve the two-fixed, center problem with the moon (fixed at L (T)) 
and initial conditions R. ,% and R. - to obtain an approximation at tune T to the 
restricted problem with initial conditions R 10 and Rj^ and moon initially at L(0). 

The analysis for a system rotating about any point other than the bary- 
center is carried out in a similar way, but the algebra is more complicated. 

The origin of the rotating system is to be the point A, defined by Eq. (4), with 
a and y determined from Eqs. (58) and (60). 

In Fig. 3, the vector A and the original and modified initial conditions are 
shown in the rotating system. 
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Again, it is seen that the two-fixed center problem, with primed initial con- 
ditions and unprimed positions of earth and moon, is related to that with unprimed 
initial conditions and primed positions of earth and moon by a rigid rotation which 
is the rotation part of the transformation carrying the rotating system into the fixed 
system. It must be remembered, however, that unlike the barycenter B, which 
may be regarded as a fixed inertial point, A is an accelerated point in inertial 
space, so that more than a rotation is required to transform back from the rotating 
system to the fixed system. In Fig. 4, the system rotating about A is shown at 
t = 0 and t = T. 

It is now easy to see that the translation required to complete the transforma- 
tion to axes moving with A, but with fixed directions, is a translation from A to A . 
Actually, this translation need not be considered further because it is desired to 
find modification in the initial conditions relative to the earth rather than relative 
to A. 


Referring again to Fig. 3, it is seen that the primed positions of the earth 
and the moon define a line parallel to that of the earth and moon at time T in the 
fixed system. Thus, just as in the barycenter case. 


and 




Ae 


= P. 


( 66 ) 


To obtain Ae, one may note that Ae is obtained by a rotation of E through CO T 
about A and that this Ae is just the negative of a rotation of A through co T about 
E. The vector A, relative to E, is given by 


A„ = A + 


n r — jtx < / 

“E ~ M+M “ 


L = (O' + 


M + M 


7 ) L + y L , 


(67) 


and the change in A_, induced by a rotation of A„ through co T about E is given 

v “hi — hi 

by 


Aa = (Oi + - 
E M + M 


7 ) (L (T) - L (0) ) + y (L (T) - L(0) ) 


(68) 


= - AE , 

so that finally, 

Rio = E 10 + + TJ7T7 ) (k (T) - L (0) ) + y (L (T) - L (0) ). (69) 

As before, P, which may now be regarded as Rio in.the fixed system, is unmod- 
ified. The two-fixed center problem, with Rio and Rio as initial conditions with 

the moon fixed at L (T) relative to the earth, should produce, at time T,.a good 
approximation to the restricted problem, with initial condition Riq and R 10 and 
the moon initially at L (0), provided T is small enough so that the second and 
higher order time derivatives of J„ produce a negligible effect. 
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PRELIMINARY NUMERICAL RESULTS 


The parameters a and 7 have been determined for a lunar orbit with the 
following initial conditions: 

x 1Q = -37163.638 km 

y. . = -56452.867 km 

J 10 

- -30844.317 km 

x Q = -0. 65536162 km/sec 

= -2. 7369109 km/ sec 

z, „ = -1. 0459904 km/sec 


The distance of the vehicle from the earth is about 11.6 earth radii, and it has 
a speed of about 3 km/sec. For these conditions, the values of a and 7 are the 
following: 

a = -6. 2611792 xlO" 4 

7 = 0.28110731 hr 


The two-fixed-center calculation with the initial conditions modified for 
evaluation of the position and velocity of the vehicle at 23, 33, and 53 hours was 
compared with the base orbit at 23, 33, and 53 hours respectively. The devi- 
ations in position of the two-fixed-center calculation from the base case are shown 
in the table below. Included in the same table are the deviations of the corres- 
ponding Kepler problem from the base case. 


Time 

Dist. from 
Earth 

Deviation 

Two-Fixed 

-Center 

Kepler 

23 hr 

35.3 ER 

A x 

Ay 
A z 

144 km 
132 km 
33 km 

17 0 km 
200 km 
10 km 

33 hr 

42.1 ER 

A x 
Ay 
A z 

262 km 
155 km 
142 km 

430 km 
250 km 
30 km 

53 hr 

52.7 ER 

Ax 
Ay 
A z 

1300 km 
1080 km 
993 km 

1970 km 
1100 km 
110 km 
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It can be seen from the table that the deviations resulting from the use of the 
two-fixed -center problem are slightly smaller than those of the Kepler problem. 

It is desirable to obtain much smaller deviations than these, but, because a and 
y used are determined only from the initial conditions, one could hardly expect 
better results. The use of one of the more sophisticated methods for determining 
or and 7, outlined earlier, should lead to considerable improvement. As noted 
earlier, these methods would render <y and y dependent on time as well as on the 
initial conditions. The smallness of the deviations (all are under \%) indicates 
that times of at least to 60 hours could be used without prejudicing the validity of 
the approximation. 
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SUMMARY 


f (- 

A recursive process for generating multi -variable ortho- 
gonal polynomials is developed. 


I . INTRODUCTION 

Let [P 0 ,XO' )], [P 1 ,X(P 1 )], .... t0 n ,X(0 n )] be a collection 
of tabular points, where 3 = (t Q ,t-^, . . . ,t m ) and 3^ = (t^Q,t^^, . . . ,tj m ) . 
That is. 
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» 0 " ^ t 00» t 01* * * * ^Om) 
e i = ^ t 10» t ll»*** ,t lra 


where the first subscript denotes the particular point and the 
second subscript denotes the particular variable. The least 
squares problem for several variables is that of finding a 
polynomial 

N 

W 0) + A 'i Cp l (e) + ••• + A N cp N (0) = Z 0 A j Cp j (3) 


such that 



N 

Y A cp 


(3 )] 


2 


is a minimum. Sufficient conditions for the existence of this 
polynomial were developed in Cl]. 


Now suppose we let 


N 


p <V*i V - - jS o Vj < 0 i n • 

A necessary condition that this be a minimum is that — dJL- = 

s A o 


9 F 




F =0. This yields the system of equations 


N 


‘oVo + A i ,p i' p o + • • • + 'Wo " x * cp o 

Vo 5 i * A l^l + • • • + Vft ’ x **i 


A.cp.cp-T + 
0 ON 


Vi^N 


+ = X,tP N 
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where X = [X( 0^) ,X( 0^) , . . . ,X( 0^) ] and cp^ = [cp^ ( 0^ ) ,cp^ ( 0^) , . . . ,cp^ ( 0^) ] , 

j = 0,1,. ..,N. Obviously, if cp (0), cp (0), ..., qp (0) are chosen 

0 1 N 

so that cp«cp. = 6 . ., then the problem is greatly simplified. 

I J * J 

II. The Recursion Process 


Let g , g , ...» g be the following set of vectors: 
0 1 m 

g “ ( t , t ,...,t ) 

6 0 00’ 10 no 

g ! = ^01^11* ’ * * #t nl* 


g — ( t ~ , t t • • • f t ), 

m 0m lm nm 


and define the vectors g * , g', ...» g* as follows: 

0 1 m 


g Y = ®Y ' (g Y ’ e 0 )e 0 * <g Y* e i >e i 'S’Vl’Vl’ 

Y = 0,1, where e = g’/llgjand (g ,e ) = g«e , k = 0,l,...,m-l. 

y Y 11 Y" v k Y K 


Then the vectors e^ , e^, ..., e^ form an orthonormal collection. 
Notice that 

- i/||iJ|Ci k - <i k .* 0 )e 0 * ‘V^i ' •" ' 'VVi’Vi 1, 

Theorem : If A^(k) = (g Y »e k )/j|g^||, then 


A (k) = [A ( -1 ) A (-1) (g ,g, ) - A ( 0 ) A (0) - A (1) A (l) - 
Y Y k Y k Y k Y k 


... - A (k-l)A (k-1)], for k = 0,1,..., Y -1. 
Y K 
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Proof: iMk) = < S y . ® k > /|| g;|| = l/||g| g y . S U /||S^| ~ <«k- e o ) e 0 / lt«kji " " 


•” - (s k’ e k-l )e k-l / l 6 ki 1 


= <g Y ’ S k )/C li g yll ll s k|| ] ' (s k’ e o )(g Y’ e O )/[ ll g kl||| g y|| ] ‘ " 


" (g k ,e k-l’ (g y’ e k-l >,/[ ll s k|||| g yl 


A ( -1 ) A (-l)(g ,g ) - A (0) A (0) - ... - A (k-l)A (k-1) 

' J 1 ' v K Y n 


V k 


Thus , 

"e = A (-l)g - A (0)e - ... - A ( y-l) e . 

Y Y y Y 0 Y Y-l 


This theorem allows us to construct the following triangular 
array of coefficients that will be needed in later calculations; 


V' 1 ’ 




A^-l) 

A ( 0 ) 



y-i> 

V 0) 

V 11 


V - 11 

v°> 

a 3 (D 

A (2) 
3 


Notice that only the elements in the first column require any 
new calculations, since all other elements in the array can be 
written recursively using these elements and the previous theorem. 
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To show how these coefficients are to be utilized, let us 
define f (3) as follows: 

Y 

f (3) = A (-l)t - A ( 0 ) f ( 3 ) - A (l)f (3) - ... - A (y-l)f ,(3), 

Y YYY° Y 1 Y Y-l 

for y = 1,2,. ..,m, and 

V 8) - A o ( - lH o- 

Notice that each f (3) is a linear combination of the m + 1 indepen- 

Y 

dent variables t , t , . . . , t . 

0 1 m 

Theorem: If f = [f (3 ) .f. (3. ) . . . . ,f (3 )], y = 0,1, then 

Y y 0 Y 1 Y n 

f = e . 

Y Y 

? 0 = Cf O (P O ),f O ( 0 1 ), * , * ,f O ( 0 n )] 

= l> o <-l)t 00 ,A o (-l)t 10 A 0 ( ” l)t n0 :i 

“ A o (-l)[t 00 ,t 1() ,...,t n0 ] = A 0 (-l)i 0 = e Q . 

Now assume that f^ ^ = e^ Then 


f = [f (3 ) ,f, (3 ) , . .. ,f (3 ) ] 
k k 0 k 1 k n 


= [A, (-l)t - A (0)f (3 ) - A (l)f (ft ) - ... - A (k-l)f ( 3. ) J 

k Ok k 0 0 k 1 0 k k-1 0 

A (-l)t - A (0) f (3 ) - A (l)f (3 ) - ... - A (k-l)f (3 .), 

k lk k 0 1 k 1 1 k k-1 1 


A. (-l)t - A (0)f (3 ) - A (l)f. (3 ) - - A (k-l)f (3 )] 

k nkkOn kin k k-1 n 
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= y-l)Ct 0 k ,t ik ,...,t nk ] - V 0 )Cf 0 (S 0 ,lf O <B l > ’"-’ f O (S n n 

- A (i)[f (0 ),f (9 ) # ...,f.O n n 

k 10 11 1 n 

- A (k-l)[f (0 ) ,f (0 

k k-1 0 k-1 1 k-i n 

(-i)i k - A k (0)f 0 - 

(-l)i k - A k (0)e o - y l)e a = V 


= A 


= A 


Thus, if we take f (0) = cp (0), j 

j 3 

equations, then the solution is 

A = X*f 
0 0 

A = X«f 
1 1 


= 0,1,..., N, in the normal 


A 

N 



where we must have N = m. Therefore, the approximating function is 

At + A t +...+At. 

0 0 11 m m 

Now suppose we desire an approximating polynomial containing 
a constant term as well as all possible second degree terms. Ihen 
we will need to include the following vectors in this treatments 


v 2 

g m+3 


( 1 , 1 ,.. ., 1 ) 

2 2 2 
(t 00 ,t 10 ,, * #,t n 0 ) 

<t 0l ,t ll»**** t nl ) * etC ‘* 


denote typical elements of these vectors by the variables 

t t , etc., and proceed as before in deriving the A 

m+2 m+3 ' 

coefficients and the function f y ( 0 ) - 
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Summary 


A procedure using quadrature methods and combinatorial topology is des- 
cribed for computing values of integrals in n-dimensions . This offers one 
way of solving the problem of data point selection for the generation of a 
least squares approximation of a multivariate function by a linear combination 
of polynomials which are or Tnonormal over a region. 


SECTION I . INTRODUCTION 


Progress Report No. 2 on Studies in the Fields of Space Flight and Guid- 
ance Theory contained the first part of this investigation of the approxima- 
tion of functions X and t r of many variables using a least squares criterion. 
An iterative method of generating a family of orthonormal functions to was 
described. This method is now part of a computer program for approximating 
a function given in the form of tabulated data. Generalized Fourier coeffi- 
cients c± are formed using the definition 
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c ± z = Zj w k q ik^k* ^ 

k-1 

where (q.^, X) is the inner product of the function qp with the control func- 
tion X, \ is a weight applied to the k th data point. The two guidance 
functions % and T R were then approximated 

(2) 

t e * fovV®) (3) 

J = 

The methods of selection of the weights w^ and the n data points were not 
discussed extensively. This aspect of the solution will now be presented. 

SECTION II. THE GENERALIZED FOURIER COEFFICIENTS, 
ORTHONORMALIZATION, AND MULTIPLE INTEGRALS. 


Both the orthonormalization of a given set of basis functions jb^jto form 
the set jq.Q or and the formation of the generalized Fourier coefficients 

c i depend on the definition of the inner product (bp,bj) of two functions. 


Usually there are uncountably many minimum fuel trajectories which ful- 
fill a given mission if the initial conditions lie within some closed bounded 
region. These trajectories form a region R over which the functions 


T„ (t, x, i, y, y, |) and %(t, §, x, X, y, y, =) are defined. This re- 
gion R has hounds imposed by the physical aspects of the problem or by re- 
strictions on the initial values. A true least squares approximation of X 
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requires that 




Xqj_d-xdxd.yd.yd (^) d (^) dt 


Inner products are similarly defined. The accuracy of the approxi- 

mation depends on the accuracy with which the numerical values of the multi- 
ple integrals are determined* At first consideration, the computation of 
these multiple integrals seems a problem of at least the same magnitude as 
the original one of approximating a multivariate function. This is easily 
seen since 



. x r ) dx-^dxp 


.dx r - 


k 

E W 1 f ( x li> x 2i * * - x ri^ ■ E ( f ) 


(M 


If the error E(f) is zero for polynomials of degree d or less in the r 
variables, then the quadrature formula is said to be of degree d. In this 
case, the direct way to obtain the set of weights, w^, and points (x^, * 21 * 
. . . x r j_) would be to solve the non-linear, non-homogeneous, algebraic 
equations obtained from (4) by substitution of a monomial for f . 


E- i T(x Ji ) d J = f. . . f if ( x j) dJdx j 

i=l j=l J J R j=l J 


for all sets of d. such that XU d • ^d. This would involve solving a system 
(d*r) ! 

of such equations. One solution would result in a quadrature formula. 

d!r ! 
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This direct method, unfortunately, leads to more complicated problems than 
the original one of approximating inner products. However, we can now use both, 
classical and "modern 11 developments of mathematics to provide alternate methods 
of evaluating these inner products. With only a sample of trajectories, the 
problem of point and weight selection for equations 1-3 can be reduced by 
special methods to finding quadrature formulae for the simplest geometric 
figures (simplexes) in a finite dimensional space. 

SECTION III. SIMPLEXES, MULTIPLE INTEGRALS, AND QUADRATURE METHODS. 

A set of points p Q , pq, . . P r in r-dimensional Euclidean space E r is 
said to be linearly independent if the set of vectors (or elements) (pq-p 0 )> 
(P2"Po) ’ • • •> (P r -P 0 ) are linearly independent; i.e., if 

(pq-Po) + ^(^"Po) + • • • ♦ *r (Pr _ Po) * 0 

implies </ ± = <* 2 = • * • = <* r * 0 where the are real numbers. If v Q , Vq, 

. . , y v r are independent points, then the set of points p of the form 

p* “ ^o v o + *1 Vq ♦ • • • 4 ^r v r 


where 



i s 0 


and 


i ■ 0, 1, . . .> r 

is called a simplex with vertices vq. For a given point p in the simplex, the 
o(q are called the barycentric coordinates of p. 
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Any simplex which has S (S^r)of the v^ as vertices is a proper face of 
the original r-dimensional simplex. Two r-dimensional simplexes are properly 
situated if their intersection is a common proper face or the null set 0. 

A finite set of properly situated simplexes is called a complex . The set 
of minimum fuel trajectories representing possible disturbances in the state 
of the vehicle and still completing the mission are given as tabulated data. 
Along each trajectory, the values of the parameters are given as 7-tuples 
with the associated values of the two control functions and T^. These 
7- tuples define the region R for the purpose of a least squares approximation 
using numerical methods. The complete set of tabulated trajectories defines 
a complex. This complex C is an approximation of the region R. 


By decomposing the complex C into properly situated simplexes, integrat- 
ing over each simplex, and finally summing the values of the multiple inte- 
gral over each simplex, an approximation of the integral over the region R 
is obtained. In our particular case, we wish to find weights wi and points 
Pi = (ti, g, x i , x-p y ± , y i , | ± ) such that 



= Z7 w iP d (ti> W 1 '*!' FI* yi' mi) 


i-1 


where is a 7-dimensional simplex and is a polynomial of degree less than 
or equal to d in the seven variables . 


Integration over each simplex can be an iterative task for a computer 
by the use of one quadrature formula of a given degree d. The number of 
points at which a function must be evaluated for use in a (2m-l) order quadrature 
formula that is valid over an r-simplex is equal to m r in most cases. For 



195 


example, with the immediate problem at hand, if we wished to integrate a 5 th 
degree polynomial in 7 variables exactly over a 7-simplex, 3 T (=2l87) tabu- 
lated points would be needed for each simplex. If the complex C were com- 
posed of only a few simplexes, the use of such elaborate formulae could possi- 
bly be justified in terms of computer time taken and the accuracy of the re- 
sults obtained. However, the parameters and the control functions are rea- 
sonably smooth, indicating that simpler quadrature formulae requiring fewer 
data points may be used. This would be especially true if the euclidean dis- 
tances between the vertices of the simplex S r are small and the values of 
the parameters do not change rapidly within the simplex. 

Recent work by Stroud (Ref. 6) would indicate that it is possible to 
find formulas requiring far fewer than m r points for (2m-l) -degree integra- 
tion. A third degree formula for a r-dimensional simplex was developed us- 
ing only 2r+3 points and not 2 r points. Unfortunately, at this time, there 
seems to be no general theory for the generation of these simpler quadrature 
methods . 

The hypervolume A r of a simplex S r with vertices v ± . (x ±1 , x i2 , ■ • • x ir ) 
is required in the development of quadrature formulas. This hypervolume is 
easily computed in the form of the absolute value of a determinant 
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For an extensive explanation of this formula see Ref. -.5* This same determi- 
nant may be used to test the independence of the vertices (points) v^, 
i ■ 0, 1, . . . r. 

Quadrature formulas may be developed to give the exact value of the in- 
tegral of a degree d polynomial in r variables over a r-simplex S r . (See Ref. 

1, 2, 3, 6, 7, and 8.) However, for computer use, an affinely symmetric 

* 

formula is desirable. In this type of formula, the weight w^ for the point 
Pj. does not change when the r-dimensional space containing the simplex is 
affinely transformed. In other words, if w^ is the weight associated with 
the point pj_ in the simple^ S r , then is associated with the point 


a p ± . y 

(6) 

T(Pi) 

(7) 


* 

where the points p^ and p^ are written as column vectors, A is a non-singular 
matrix of real coefficients, and X is a column of constants. Equation ( 7 ) is 
(6) written with the affine operator T. is in the simplex «T(S r ), the 
set of transformed points of S f . An additional requirement may be imposed 
on a quadrature formula for a simplex. The points p^ used in the formula all 
must lie within the given simplex. This restriction is justified for two 
reasons: 

(1) The function may not exist outside the simplex. 

(2) Since any point within a simplex is determined by its barycentric 
coordinates, a simple computer routine can be used to find the quad- 
rature points from the vertices of the simplex. 


An affinely symmetric formula of third degree can now be given. Let the 
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r+1 vertices of the r-simplex be v Q , v-^, . . . , v^,. The barycenter B 

of S r is defined by 

b ■ ^ t>r 

The -hypervolume of S r isAr. Then 

f fdv . = *r E f (u^ ) ♦ c r f(B) 

J S r r o 

where 

p - (r+3) 2 A 

^ - 4(r*l)(r*2) 

c - -( r - > - 1 ) 2 A 

r t(r*2) Ar 

u ± - - 2 - v. + Ili B i - 0, 1, . . ., r 

1 r + 3 1 r+3 

The formula is exact when f is a 3rd degree polynomial in r variables. 

This means that the values of inner products of the types (x-j™, x^ 11 ), 

(m»n = 3; i, j = 1, . . . r) will be exact over each simplex S r in the region 

R. If the i nn er product is of the type (x d , f) ; i s 0, 1, 2, where f is not 
a polynomial of degree 3 or less, then there will be an error due to the 
quadrature formula. There has been little error analysis available for 
quadratures involving functions of many variables. However, the errors 
arising by using such a definite procedure as the above are usually much 
less than if a simple sum of products had been used to approximate an inner 
product . 

The approximation by the use of quadrature formulas is an approximation 
over an entire region R and not over a finite set of points as in a least squares 
method such as normal equations or orthonormalization of vectors. 
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SECTION IV. RECOMMENDATIONS 

In order to arrive at a suitable general algorithm for the approximation 
of control functions for Saturn class vehicles, it is recommended that the 
studies contained in this report and Part I be continued. This continuation 
should include the following particular areas of effort: 

1. The development of a method, suitable for computers, for finding the 
vertices of all the properly situated simplexes in the region defined 
by the tabulated minimum fuel trajectories. 

2. The implementation of available information from the calculus of vari- 
ations and multivariate functions to determine the boundary of the 
region over which minimum fuel trajectories are defined for a parti- 
cular mission. 

3. The comparison of the accuracy of an approximation using quadrature 
methods to define the inner product of two functions with the usual 
method using sums of products of the values of the two functions at 
arbitrary points. 

4. The study of possible methods of directly producing a rational approx- 
imation from a polynomial approximation or a partial sum of a series. 

5. The investigation of direct substitution of a polynomial with undeter- 
mined coefficients for the control functions into the Euler-Lagrange 
equations; the goal being the determination of the coefficients which 
will minimize the fuel consumption for a particular mission. A gener- 
alization of the two point boundary problem would be needed with in- 
equality constraints or "initial conditions" satisfying some inequality. 

6. The exploitation of analog computer methods, which may possess advan- 
tages in terms of time and money, in the areas of both Chybechev and 
least squares approximations deserves renewed effort. 
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SUMMARY 


/ y *it 


The purpose of this work is to obtain simple, formal functions 
which approximate, in some sense, the steering and cutoff functions 
derived in the Adaptive Guidance Mode. The approach taken in this re- 
port is to use linear programming techniques to fit linear combinations 
of known functions or ratios of such functions to a set of tabulated 
values of the steering and cutoff functions* 


I. INTRODUCTION 


This report describes the use of linear programming techniques 
to approximate the steering and cutoff functions for the implementation 
of the Adaptive Guidance Mode. [5] [8] [10] This approach to the approx- 
imation of the guidance functions is basically a multivariate curve- 
fitting problem. Values of the steering and cutoff functions are 
tab ula ted for a representative set of points on minimum fuel trajectories 
and then, formal functions are sought which approximate these tabulated 
functions according to some criterion. When this criterion is Lj (min- 
imized sum of absolute deviations) or L» (Chebyshev or minimized maximum 
deviation), linear programming may be used to determine the approximating 
functions. An analysis has been made of the case in which the approxi- 
mating functions are polynomials. Studies have been initiated on the 
use of ratios of polynomials for the appr oxim ating functions. 
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Section II of the report contains scree results in the theory of 
linear programming which are included as background for the later dis- 
cussions • 

Section III contains a statement of the curve-fitting problem. 

In section IV, the case of the Lj approximation of the guidance 
functions by a polynomial is considered. This case is included for com- 
pleteness and for purposes of comparison with the L» case, which is 
of more interest in practical applications. 

The Ip® approximation of the guidance functions by a polynomial 
is considered in section V. The problem is also formulated so that a 
linear programming routine can be used to find that polynomial, if it 
exists, which approximates the tabulated function to within predeter- 
mined tolerances at each data point. 

Section VI contains a discussion of peculiarities of the curve - 
fitting problem which cause slow convergence of the linear programming 
method. Recommendations for improving this speed of convergence are 
included. 

Numerical examples of the Lj and L» approximation of the steer- 
ing function by polynomials are given in section VII. 

Section VIII contains a brief discussion of experiments done using 
alternative methods for choosing the pivotal elements in the simplex 
algorithm for linear programming. The purpose of this work was a fur- 
ther increase in the speed of convergence for the simplex method. 


II. THEORETICAL BACKGROUND 


Linear programming problems which arise in curve-fitting can often 
be solved more readily in their dual form than in the original primal 
form. The basic properties of dual linear programming are therefore 
summarized in this section. 

The Duality Theorem states that if the primal (dual) problem has 
a finite optimum solution, then the dual (primal) problem has a finite 
optimum solution and the extrema of the respective objective functions 
are equal. If the primal (dual) problem has an unbounded optimum solu- 
tion, then the dual (primal) problem has no feasible solutions. 

If a bounded optimum solution for the primal problem exists, then 
the solution of the dual problem can be obtained by solving the primal 
problem, and vice versa. The desired solution can therefore best be 
obtained by solving the simplest of the primal and dual problems. 
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There are several pairs of dual linear programs. [1] [2] [3] The 
most familiar pair consists of the ” canonical" or "symmetric” dual pro- 
grams. Goldman and Tucker [3] point out that the other problem pairs 
are essentially no more general than this canonical one. Several pairs 
of dual linear programs are discussed below. In practice, the appro- 
priate pair must be chosen to fit the special requirements of the par- 
ticular curve -fitting problem. 


(A) The canonical pair of dual linear programs is stated as follows: 


Primal problem (dual problem): 


( 1 ) 


Minimize the objective function, f a CX , subject to the 
constraints, AX i b and X * 0 , where X is an n -component 
J column vector of unknowns, C is an n-component row vector, 

J A is an m X n matrix, and b is an m -component column 
^vector. 


Dual problem (primal problem) : 

(Maximize the objective function, g » Wb , subject to the con- 
Jstraints, WA < C and W s 0 , where W is an m-component 
(row vector of unknowns and A , C , and b are as defined 
[in (1). 


To solve (1), the constraints, AX > b , are converted to equalities 
of the form, (A , - 1)^X^ " b , where X is an m-component column 

vector of non-negative "slack” variables and I is the m x m iden- 
tity matrix. The linear programming problem is then stated as, 

(Minimize the objective function, f * CX , subject to the 
J constraints , 

- b , X * 0 , X*0. 

To show that the solution of (2) can be obtained by solving 
(3), we let Xq and ^ denote the optimum solutions of the primal 
and dual problems, respectively. Let B denote the optimum basis 
of (A , - I). Xq is an m-component, column vector of basis var- 
iables (its elements correspond to those columns of (A , - I) which 
occur in the basis B ). Now let % denote an m-component, row 
vector, each element of which is an element of C corresponding 
to a basis variable. Consider the vector W 0 ■ C^B" 1 • B" 1 A is an 

m x n matrix, each element of which is an element in the simplex 
tableau for the solution of (3), and Xq is the optimum solution 




of (3); therefore all the shadow prices are non-negative* Hence, 

C - CqET 1 ! * 0 and 0 - C^ET 1 («l) i 0 , Therefore W 0 A £ C and 
Vq > 0. Hence is a feasible solution of (2)* Furthermore, 

BXq ■ b j therefore W 0 satisfies C^Xq » • It then follows 

from the Duality Theorem that Wq is the optimum feasible solution 
of (2). 


The unsymmetric pair of dual programming problems can be 
stated as follows: 


Primal problem (dual problem): 


(U) 


Minimize the objective function, f ■ Cl , subject to the 
constraints AX « b and X i 0 • 


Dual problem (primal problem): 


(5) 


Maximize the objective function, 
constraints, WA £ C , where W 


g ■ Wb , subject to the 
is unrestricted in sign. 


This pair of problems can be obtained from (l) and (2) by expressing 
the equality constraint AX ■ b as a pair of inequality constraints, 
AX i b and - AX i -b * By a proof similar to that for the canon- 
ical pair of dual problems it can be shown that the solution CqB - 1 
of (5) can be obtained by solving (U), where C*> and B are as 
previously defined. 


The pair of dual problems for a linear programming problem 
with mixed constraints can be stated as follows: 


Primal problem (dual problem): 


Maximize the objective function, 
(6) < t° the constraints. 


JA ll X 1 + A 18 X S i bj , 
(unrestricted in sign. 


^si^i + ^sa^s 


GjXj + CgXg , subject 
hg , Xj 2 0 and Xg 


Dual problem (primal problem): 


(7) 


I Minimize the objective function, 
Jto the constraints. 


g 


IWjAjj + WgAgj S Cj , 

l unrestricted in sign. 


+ ^a^as 


Wjhj + Wgfc^ , subject 
Qj , Wj > 0 and Wg 


The solution, Gq ET 1 , of (7) can again be obtained by solving (6) 
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in. THE CURVE-FITTING PROBLEM 


The purpose of this work is to obtain simple, formal functions 
which approximate, in some sense, the steering and cutoff functions for 
a missile on a minimum fuel trajectory. The approach taken in this 
report is to use linear programming techniques to fit linear combina- 
tions of known functions or ratios of such functions to a set of tabu- 
lated values of the steering and cutoff functions. The studies were 
performed for a flat, two-dimensional earth with no atmosphere, a point 
missile, and constant fuel flow. The steering function, h , and the 
cutoff function, T , are functions of the six independent variables, 
x and y (rectangular space— fixed position coordinates ) , x and y 
(velocity coordinates), F/m (thrust acceleration), and t (time). 

The approximation of the steering and cutoff functions is studied 
by considering the general problem of approximating a function, 
f(z) , whose value is known at n points, Zj , » in an m- 

dimensional space, by a function, P(z) * of a given form in the com- 
ponents of the vector z • When the criterion of "best fit" of P(z) 
to f(z) is Lj , i.e., the sum of the absolute deviations, 

n 

X jPCZk ) - f (.%,)] , is minimized, or If» , i.e., the maximum abso- 

k=l 

lute deviation, max |P(^) - f (z^ ) I » i s minimized, then the 
k=l,ljn 

function P(?) can be determined by linear programming techniques. 

The function P(z) is assumed to be of the form, 
m m 

P(Zk) ”A 0 +ZA t Pj(z kl )+ XA tJ P 1 j(z kl ,z kj ) 
i»l i>j«l 

m 

+ Z|ti ’ * k J ) * 


•••) 
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where the Aj 's , A {] <8 , ••• are unknown coefficients to be deter- 

mined, ^ t is the i-th component of ~z at the k-th data point, and 
Pj (Zfc . ) , axB predetermined functions of z, , 

(zj z i ) , ..., respectively. In sections IV-VII, P(z) is a polynom- 
ial, i.e., Pt (*ki ) “ j (*k i > j ) “ *k i Zk j > etc. 

Various pairs of dual linear programming problems are formulated 
in sections 17 and V for the curve-fitting problem. The number of 
constraints (other than those requiring variables to be non-negative) 
in one problem of the pair is usually of the order of n , the number 
of data points. The number of constraints for the other problem of 
the pair is of the order of the number of unknown coefficients in 
P(z) . In practice, n is very much larger than the number of co- 
efficients. The time required to compute the optimum solution is 
approximately proportional to the cube of the number of constraints. 
Hence, a significant decrease in computation time can be obtained by 
solving that problem of the pair of dual problems in which the number 
of constraints is a function of the number of unknown coefficients in 
P(S) . 


IV. THE Lj APPROXIMATION 


In an 1^ approximation, the function, P(z) , which minimizes 

n 

the sum of absolute deviations, X | P(z, ) - f(z^)| » is sought. For 

k=l 

m 

simplicity, only functions of the form P(5^ ) ■ Aq + l A |2*i "in 
be considered. i=*l 

The curve -fitting problem can be restated as the following linear 
programming problem. [11] 
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Minimize the objective function, 
constraints. 


n 


X (€ k + 6* ) 

k=l 


» 


subject to the 


( 8 ) < 




+ Y,k i z ki - € k + 6 k - f(Zfc) t 


i-1 

€ k * 0 and $ k * 0 , 

where the A t 's are unrestricted in sign. 


(k - 1 , n) 


In order to solve (8) by the simplex method or the revised simplex 
method, the unknown coefficients, A t , must be expressed as the dif 
ference of two non-negative unknowns, i.e., A. ■ a t - b t , where 
a t at 0 , b t it 0 . With this substitution, (8) becomes: 


n 


-ion, Z (€ k + 6 k ) t 


Minimize the objective function, L> (€ k + 6 k J , subject to 
the constraints , k**l 


m 


" f K ) * 


(9) “S Sq - b^ + X (a, - b, )z kt - € k + 6 k 
i-1 

€ k ^ 0 , iO 

and a, > 0 , b t * 0 for i*0,..,m 


(k - 1 , .., n) 


The computation time required to obtain the optimum solution of (9) is 
approximately proportional to n* , the cube of the number of constraints. 
Since n will be very large in practice (i.e., of the order of 3000), 
the computation time will be lengthy and can be reduced by solving the 
dual problem to (8). 


By rewriting (8) in matrix farm and applying (6) and (7), the 
following pair of dual problems is obtained. 


Primal problem : n 

Minimize the objective function, E - ^ (€ k +8 * ) * subject to 
the constraints . k-1 
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1 • • * ®x ,m 0 . # • 0 1 0 • • 0 


• • • 


jjL *jji ®n2 • • • z n,m 0 0 • • •-! 0 0 * . lj 


where ^ ^ 0 ^ for k * X | *t| Hj 

and Aq , \ sure unrestricted in sign. 

Dual problem : n 

Maximize the objective function E' « l f (4K > 
the constraints. 


^0 


f(0 

K 


f(%) 

• 

m 

• 

# 



f( z n) 


— j 

*n 




*n 



k=l 


subject to 


1 1 

z n 

*12 


• • » 


1 

• • z nl 

* * z n2 


L z i,. % 


t* • • • z n,mj 


% 


0 



0 

• 


• 

• 

- 

• 

_v 


_0 


and 


-1 o ... 0 


”l 


1 

0 -1 ... 0 







Ua 


1 

0 0 . . .-1 


• 

< 

• 

1 0 ... 0 


• 


• 

0 1 ... 0 







”n 


1 

_ — 

_ 0 •••*!_ 




where the Uj 's are unrestricted in sign. 


Here the dual problem has m+ 2n + 1 constraints, so that the computa- 
tion time is proportional to (m + 2n + l) 3 . Hence, the situation is 
not improved by considering the dual problem directly. However, if a 
computer program for the simplex algorithm for bounded variables is 
available, the following equivalent form of the dual problem, with 
m + 1 constraints, can be solved instead with resulting savings in 
time. [12] 
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Let v t * u, + 1 for i - 1 , 2 , 


• • •> n « 

n 


k=l 


Maximize the objective function, E' ' - I { f )v k - f (4 )i , 

n 
n 

i-. 


subject to the constraints, 

1 X • • • 1 


*11 231 


J nl 


*!,■ *3,» • • • z n,m 


’n 


k=l 
n 

k=l 
n 

y z *. 

« j 


where 2 2 v t * 0 for i - 1 , 


, n * 


V. THE Ip° APPROXIMATION 

In an Ip approximation, the function, P(z) , which minimizes 
the maximum deviation, max |P(z k ) - f(zjc)l , is sought. 

k*l, • ,n 

By introducing a positive number € , the linear programming prob- 
lem can be stated as : 

Minimize € subject to the constraints, 

P(^ ) - f (z* ) < € 

P(4 ) - f(^ ) 

m 

Assuming P(zJj ) • Aq + I A t a* j , the problem becomes : 

i=l 

Minimize € (or maximize - €) subject to the constraints. 


for k * 1 , • n 


m 


(10) 




+ ^ Aj z k j - 


€ « f & ) 


for k = 1 , •«, n and 


i“l 
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Bt 

- Aq - I A lZkl - e<-f (z*) 

(10) i«l 

where the A t *s are unrestricted in sign and € is non-negative. 

This problem has 2n constraints. Before solving (10) by the simplex 
method, the unknown coefficients A t must be expressed as the differ- 
ence of two non-negative unknowns, i.e.. A, • a t - b t . As in the 1^ 
case, the computation time may be decreased by solving the dual problem. 
By considering (6) and (7), it can be seen that the dual problem for 
(10), stated in matrix form, is the following problem with (m + 2) con- 
straints, In deriving the dual problem, all the variables in (10) are 
considered to be unrestricted in sign, since the form of the constraints 
in (10) ensures that € is non-negative. 


n 


r 


y 

>ion, L, 


Minimize the objective function, ^ [f(zjf )^ - f (zjj )v k ] , subject 
to the constraints, k=l 


01 ) { 


1 

1 


. . 1 

. • 1 


1 

-1 


. . 1 

. . -1 


& nl "*ii 


-z 


nl 


13 »■ * * z n,ra “*!,«•• “®n,m_ 


"n 


v i 


1 

0 


3n_ 


\and ^ 0 j i 0^ for 1c * X ^ • • • j n • 


Similar techniques can be used to determine the polynomial, P(z) 
of a given form which approximates f(z) while keeping the deviation, 
|p(zit) “ f(z k )| , at the k-th point within a predetermined tolerance 


€ k , for k * 1 , ,,, n . Assuming that 
introducing \ = l/€ k , for k ■ 1 , ,,, 


m 


“ Ao 


. I 


i=l 


3^1 


and 


n , the problem becomes that 
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m 

of determining Aq , .., A^ so that max \A 0 + x*,*. - f(£)l 

k 5 *! , * »,n i*l 

is minim ized. If this minimum value does not exceed 1, then the devia- 
tion of P(z) from f(z) is within the desired tolerance at each data 

m 

point. Defining Aq + L A 1 z* , - f (z* ) - \ ^ » whejre d„ i 0 , 

i-1 

i 0 , and introducing the non-negative variable t , the problem can 
then be expressed as the following linear programming problem. 


( 12 ) 


Minimize the objective function, ^ 
the constraints, 

m 

dk - e* - \k 0 - \ L k t t kt - \ f(4 ) * 

J i«l 

> d 1 +e 1 +t-d k - e* a 0 
dk a 0 , e* a 0 


+ t , subject to 


(k « 1 , n) 


a 0 , and A<, , A^ are unrestricted in sign. 


If the objective function, + e 1 + t , for the optimum solution 
exceeds 1, then there is no approximating function of the given form 
satisfying the given tolerances. In solving (12), the coefficients Aj 
should be expressed as A t ' • a t - b t , where a { a 0 , bj i 0 . The 
problem has (2n-l) constraints. There appears to be no saving in com- 
putation time in solving the dual problem* 


An alternative formulation of the problem follows: 


Maximize y subject to the conditions, 

|p(z* ) - f(^ )| * € k - y , 

i*e. 

Y + P (**)<€*+ f(^) 

and (k ■! , n) 

y ~ PC^k ) ^ ^k ) 

If the maximum value of y is non-negative, then a polynomial of the 
desired form satisfying the given tolerances has been found* Assuming 
m 

P(Zfc ) * + t * 

i**l 


the problem then becomes: 
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Maximize y subject to the constraints 


m 


y 

+ Aq + ZjA 1 z kl £ € k +f(z k 


) and 


( 13 ) i 


i=l 

IQ 

-I 


(k * 1 , .., n) 


V - Ao -Lk i z kt * € k - f(%) 
i»l 


where y , Aq , . ., \ are unrestricted in sign. 

The dual problem to (13) has only (m + 2) constraints, and hence should 
be solved in preference to (13). 
n 


/" 


Minimize 

k-1 

constraints. 


Llie* + 


(Hi) 


f(Zfc )]u k + [€ k - f(^)]v k } subject to the 

1 
0 

0 


J— — 


- *“ 

1 • ♦ 1 1 • • 1 


% 

1 • . 1 -1 • . -1 


% 

• 

• • z nl “®ii • * ^nl 

• • • • 

* • • • 


• 

u n 

v i 

• 

• 

n • * z m -z. • • —z 

L 9 n > m 1 Il,m_ 


_ V n_ 

id Ujf ^ 0 , v k £ 0 , for k « 

i , 

•*9 


VI. SPEED OF CONVERGENCE 


The special structure of the right-hand vector in the constraint 
equations in (11) and (llj.) causes slow convergence of the simplex algo- 
rithm. In this section, the steps of the algorithm are analyzed for 
this special case, and a transformation of the original problem, which 
increases the speed of convergence of the algorithm, is discussed. The 
notation in this section is that which is frequently used in discussions 
of the simplex method and has no connection with the notation of the 
proceeding sections. 

Let the following array denote the simplex tableau, where the 
dk 's are shadow prices ((^ - - C* ) , the v, 's are activity levels. 
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and z is the value of the objective function. 


z 

di 

• • • d 2n 

V 1 

Xll 

• • • x l,2n 

v a 


• • • x 2,2n 

• 

• 

• 

• 

• 

• 

v .+ i 

,i 

• • • x m+l,2n 


If the k-th column of the tableau is introduced into the basis, then 
the column in the tableau corresponding to the j-th column in the basis 
is eliminated from the basis, where j is given by 


7 ) / X J* 


min (vj /x lk ) 

i»l,..,m+l 




> 0 


The change in the objective function is - Od* • Since is posi- 

tive and 9 ^ 0 , the objective function is decreased only if 9 > 0 . 
Initially, the vector of activity levels is the column vector (1,0,.., 0), 
Hence the value of the objective function will be improved initially 
only if x lk £ 0 for i ■ 2 , 3 » .., m+1 and Xj k >0 . If the 
problem has a bounded optimum solution, then Xj k > 0 if x, k £ 0 lor 
i * 2 , ... m + 1 . Therefore the probability that 0 > 0 is 2 
when the vector of activity levels is the column vector (1,0,.., 0) and 
the signs of the x., »s are assumed to be random. If x lk > 0 ini- 
tially, for some i > 1 , then the new vector of activity levels is 
still (1,0,.. ,0) . However, when x lk ^0 for i-2,.., m+1, 
then the new vector of activity levels will contain at least two posi- 
tive elements. By the same argument, the probability that the objec- 
tive function will be decreased at the next iteration is at least 
2 ~mn * ^he same argument is repeated until the optimum solution is 
obtained. The expected number, N , of iterations required to reach the 
optimum solution will thus satisfy the inequality. 


2 ra £ N £ 2 m_1 + 


2 ° < 2 mhl 


in the special case where the right-hand vector is (1,0,.., 0) . This 
is much larger than the estimated 3(m+l)/2 or 2 (m+1) iterations 
usually required to reach an optimum solution# 

Before solving the problem by the simplex method, a transformation 
can be applied to the problem in the following manner in order to elim- 
inate the zero elements in the right-hand vector which slow the conver- 
gence. The original problem is to minimize CX subject to the 




214 


constraints, AX * B , and X i 0 , where b is the column vector 
(1,0, ..,0) • Let D be the (m + 1) x (m + 1) transformation matrix 

1 0 0 . . 0 

1 1 0 . . 0 

1 0 1 . . 0 

• • • • 

• • • • 

10 0. . 1 _ 

Consider the transformed set of constraints, DAX » Db . 

Each component of the right-hand vector, Db , is 1 . Since D is 
non-singular, the problem of minimizing CX subject to the constraints, 
DAX ■ Db , and X £ 0 has the same solution as the original problem. 

The solution, W 0 = GoB -1 , of the original problem, where B is the 
optimum basis, can be readily obtained from the solution W 0 ' of the 
transformed problem since (DB) -1 is the optimum basis of the trans- 
formed problem. Hence, ■ Cp (DB) -1 * Ct ) B” 1 D” 1 ■ V^D” 1 and there- 

fore W 0 -W 0 ' D . The final simplex tableaux for both problems are the 
same, since if P k represents the k-th column of the final tableau for 
the original problem is B“ 1 P k , and that for the transformed problem 
is (DB)" 1 (DP k ) . 


VII. NUMERICAL RESULTS 


This section contains the results of preliminary tests run at UNC. 
In these tests, only linear terms were included in the approximating 
polynomial for x • 

P(Zk ) - A 1 + A^ + Agy k + AgXy + + A;. (F/m^ + Agl* . 

Five evenly spaced data points on each of five trajectories were used. 
Difficulty in obtaining an optimum solution was experienced in the 
early stages of the tests because of the wide range in the magnitudes 
of the variables x , y , x , etc. However, after scaling the elements 
of the constraint matrix so that all elements were of the same order 
of magnitude, the optimum solution was successfully obtained. 

Lj approximation : The results of this calculation are included for 

comparison with the L» case. 


215 


i 

A, 

0 

-203.736 

1 

.236515 

2 

.0576878 

3 

U.75139 

1* 

-37.0782 

5 

2663.38 

6 

-.721*890 


|p(zk ) - x(4 )! 


5.1081*3 


max|p(4) -x(4)l * *971*130 . 

Ip, approximation: Both the primal formulation (10) and the dual form- 

ulation (11) were tested. For the dual problem the effect of the 
transformation discussed in section VI was tested. The values of the 
coefficients. A, , computed in the three different ways are shown in 
the following table. 



Primal 

Dual 


transformed 

untrans formed 

Ao 

100.382 

121.861 

23.918 


.203273 

.183696 

.187623 


.0033899 

-.000799U65 

-.0010881*3 


-25.1*105 

-27.1*675 

3.71590 

A* 

-56.1931 

-51.9561* 

-65.1*972 


1*511.03 

1*189.05 

501*2.98 

A* 

-.291315 

-.219835 

-.303201* 

A 7 (€) 

.7823 

.7823 

.7710 

Value of 
objective 

.7823 

.7823 (.7882*) 

.7823 (7823*) 

function. 



30 

No. of 

not available 

21 

iterations • 


.8627 

1.709989 

Max |P-x| 

.7823 


*This value was computed by Xq *» (B* 1 b) 0 • 
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There are several points to be observed: 

1) There is quite a large discrepancy between the three solutions, prob- 
ably caused by the accumulation of round-off errors in the perfor- 
mance of the single precision, floating point arithmetic used in the 
linear programming routine. 

2) The values of the objective function agree up to four decimal places 
for the three trials. This value should be the least upper bound 
for max] P(z lt ) - x (^ ) I • However, for the dual formulation of the 

k 

problem, max|P(4) - x(^k ) I > € • For the untransformed dual prob- 
k 

lera, the discrepancy is large, 

3) The number of iterations required to obtain the optimum solution of 
the transformed dual problem was smaller than the number required for 
the unt rans forme d dual problem, as expected. However, the number of 
iterations for the un transformed problem was smaller than expected. 
This was probably caused by the failure of the elements of the tab- 
leau to satisfy the random sign hypothesis in this small test. In 

a realistic case in which the number of constraints is large, the 
hypothesis of random signs will be more nearly satisfied and the 
relative decrease in computation time for the transformed problem 
should be larger. 

In order to reduce the effect of round-off errors in the calcula- 
tion of B" 1 for the dual problem, B -1 was computed iteratively by 

V&X - Bf l (2I -BBp) , 

where B^ 1 is the i-th approximation to IT 1 and I is the identity 
matrix. In this experiment, E^ 1 was the matrix obtained by the sim- 
plex method. Only one iteration was carried out. The set of coeffi- 
cients in P(s) were computed from Wq ■ C^Bf 1 • The result is given 
in the following table. The results of the previous calculation for 
the primal problem are included for comparison. 
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Primal 

Dual 



transformed 

untransformed 


100.382 

100.380 

100.37U 

Al 

.203273 

.20328U 

.203285 

A, 

.0033899 

.0033922 

.0033922 

As 

-25.U105 

-25.1016 

-25.U120 

A* 

-56.1931 

-56.191*3 

-56.195U 

A« 

U511.03 

U5H.ll* 

U511.22 

As 

-.291315 

-.291353 

-.2913U9 

A? (€) 

.7823 

.7801 

.7823 

Value of 

objective 

function 

.7823 

.7823 (.7773) 

.7823 (.7823) 

Max|P-*| 

.7823 

.7879 

.7829 


Improvement in accuracy is apparent. The three sets of coefficients 
obtained independently agree to four figures. Max) PCz* ) - x (2* ) I is 

k 

much closer to the value of the objective function. However, the fol- 
lowing unexplained points were observed: 

1) * l B is very close to I but BE^ 1 differs considerably from I • 

2) B has a greater deviation from I than does Eg 1 B • 

3) The accuracy of A^ » Ag seems to be improved, but that of € 
decreased. 

Vni. DIFFERENT CHOICES OF PIVOTAL ELEMENTS 


In the progress report, [6j , several criteria for the choice of 
pivotal elements were discussed. The Greatest Absolute Ascent Method 
and the Modified Gradient Method were tested in addition to Dantzig's 
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usual criterion. The accumulation of round-off errors appeared to be 
greater in the first two methods than in Dantzig's method and led to 
incorrect decisions. For this reason the alternative criteria do not 
appear to possess the merit claimed by Quandt and Kuhn. [9] This dif- 
ficulty might be removed by using double precision arithmetic or by 
improving the inverse of the basis by some iterative method after sev- 
eral iterations of the algorithm. 


IX. CONCLUSIONS 


Lj and Ip® approximation problems, or modification of them, can 
be formulated as linear programming problems. Comparison of this approach 
with the least squares procedure is difficult. The choice of a method 
depends upon the requirements and characteristics of the approximation 
problem. For example, if the deviation of the approximating function 
from f(z) is to be within a prescribed tolerance, then the linear pro- 
gramming approach is straightforward and effective. An expected merit 
of the use of linear programming is the accuracy or stability of the sol- 
ution. In the least square method an ill -conditioned matrix of normal 
equations often causes trouble. Round-off error nay cause inaccuracies 
in the linear programming approach. In particular, inaccuracies may 
arise in the process of obtaining the solution of the primal problem 
from that of the dual. This difficulty can be removed more easily in 
the linear programming approach than in the least squares method. The 
scaling problem is less complex in the linear programming case than in 
the least squares case. However, the computation tine for obtaining the 
approximating function is at least as long by linear programming as by 
least squares. 

The major results are summarized as follows t 

1. Lj or Ip® approximations of f(z) by a function 

m m 

P(z) - Ao + I + IajjP,,^,, *j) + . . .. 

i-1 i>j 

can be formulated as linear programing problems. 

2. For the multivariate approximation problem the dual program 
can usually be solved more quickly than the primal program. 

The solution for the primal program can easily be obtained 
from that of the dual program. 

3. The dual program has a special structure which may cause slow 
convergence to an optimal solution. Therefore, the constraint 
matrix should be transformed (an example of such a transform- 
ation is in VI) before the problem is solved. 
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k. The data of the original simplex tableau should be scaled so 
that the data are of the same order of magnitude • 

Sin gle precision arithmetic seems to be insufficient for obtain' 
ing the solution of the primal problem from that of the dual 
problem. A large error may occur in ^ - CqB " 1 even if an 
accurate solution for the dual program is available. There- 
fore the inverse matrix should be computed in double precision 
arithmetic. The reinversicm of the current basis is also 
advisable. 
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